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STUDY  OF  PHOTON  CORRELATION  TECHNIQUES 

FOR  PROCESSING  OF 
LASER  VELOCIMETER  SIGNALS 

By  William  T.  Mayo,  Jr. 

Summary 

The  objective  of  this  contract  was  to  provide  the  theory 
and  system  design  for  a  new  type  of  photon-counting  processor 
for  low-level  Dual  Scatter  laser  velocimeter  (LV)  signals  which 
would  be  capable  of  both  the  first-order  measurements  of  mean- 
flow  and  turbulence-intensity  and  also  the  second  order  time 
statistics:  cross-correlation,  auto-correlation,  and  related 
spectra. 

This  report  provides  a  general  Poisson  process  model  for 
low-level  LV  signals  and  noise  which  is  valid  from  the  photon- 
resolved  regime  all  the  way  to  the  limiting  case  of  non-station- 
ary  Gaussian  noise.  Computer  simulation  algorithms  and  higher 
order  statistical  moment  analysis  of  Poisson  processes  have 
been  derived  and  applied  to  the  analysis  of  photon  correlation 
techniques.  A  Dual  Correlate  and  Subtract  frequency  discrimina¬ 
tor  technique  is  postulated  and  analyzed.  Expectation  analysis 
indicates  that  the  objective  measurements  are  feasible.  Error 
analysis  for  the  mean-flow  case  indicates  that  practical 
transonic  wind  tunnel  measurements  are  possible  with  100-1000 
times  less  light  than  is  required  for  burst-counter  processors. 

A  system  design  for  a  new  high-speed  photon  processor  for  LV 
signals  is  provided. 
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INTRODUCTION 


The  Problem 

Classical  laser  velocimeter  (LV)  electronic  signal  pro¬ 
cessing  techniques  are  sometimes  inadequate  for  detection  of 
light  scattered  by  email  scattering  particles  which  are  required 
for  following  fluid  motions.  In  other  situations  detection  of 
larger  scattering  particles  io  difficult  due  to  limited  system 
sensitivity.  Photon  counting  techniques  offer  improved  system 
sensitivity  by  allowing  velocity  measurements  to  be  made  even 
when  there  are  insufficient  signal  photons  available  to  define 
the  classical  scattering  signal.  Such  techniques  are  thus  appli¬ 
cable  when  the  presently  used  classical  burst-counter  and  fre¬ 
quency  tracker-techniques  are  not. 

The  general  objective  of  this  contract  was  to  provide  the 
system  design  for  a  new  type  of  photon-counting  processor  for 
low-level  Dual  Scatter  LV  signals  which  would  be  capable  of  both 
the  first-order  measurements  of  mean-flow  and  turbulence-inten¬ 
sity  and  also  the  second  order  time  statistics:  cross  correla¬ 
tion,  auto-correlation,  and  related  spectra.  This  was  to  be 
accomplished  by  extending  the  preliminary  feasibility  analysis 
developed  under  a  brief  NASA  Langley  sponsored  study*  in  early 
1974.  In  addition  to  theory,  the  system  design  would  incorporate 
judgement  based  *on  experience  in  the  experimental  hardware  devel¬ 
opment  [1]  of  a  related,  but  simpler,  photon-counting  processing 
system  designed  and  constructed  for  the  U.S.A.F.  Arnold  Engi¬ 
neering  Development  Center  to  measure  mean-flow  velocities. 


* 

The  final  report  for  that  study  (Contract  NAS1-13140) 
was  informal  and  not  disseminated.  This  report  contains 
revised  versions  of  all  the  necessary  mathematics. 


Background 


Signal  modeling.-  Earlier  modeling  efforts  have  treated 
LV  signals  for  which  the  noise  could  be  considered  as  additive 
independent,  stationary,  and  Gaussian  [2, 3, 4, 5] .  This  is  the 
limiting  case  of  stationary  Poisson  shot  noise  which  occurs  for 
visible  light  photodetection  when  a  steady  light  source  such 
as  a  heterodyne  reference  beam  [2]  or  high  background  light 
level  [3]  dominates  the  signal.  In  a  recent  simulation  of  low- 
level  dual  scatter  signals,  the  accuracy  of  the  noise  model  was 
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Figure  1.  Triply  Stochastic  Nature  of  Low-Level 
LV  Signals:  Turbulence,  Bursts,  and 
Photo  electron  Pulses. 
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extended  by  treating  it  as  a  nonstationary  Gaussian  process 
whose  variance  is  proportional  to  the  incident  optical  power. 

When  we  discuss  "noise"  in  LV  signal  detection  we  are 
usually  referring  to  the  variation  of  the  electronically 
detected  signal  with  respect  to  a  scaled  version  of  the  classi¬ 
cal  optical  (power)  signal  incident  on  the  PMT.  In  a  general 
analysis  the  classical  LV  signals  are  also  random  processes  due 
to  the  random  amplitudes  and  the  arrival  times  of  the  signal 
bursts.  Mayo  [2]  and  Adrian  [7]  have  treated  these  signals  as 
a  Poisson  process  for  steady  flow,  and  Durrani  [4]  and  George 

[8]  have  treated  them  for  the  turbulent  flow  case  as  Gaussian  in 
the  limit  of  high  particle  number  density.  A  new  book  by  Snyder 

[9]  treats  generalized  Poisson  processes  in  great  detail. 

Papoulis  [10]  provides  a  good  introductory  treatment.  Snyder 
treats  "doubly  stochastic  Poisson  processes."  These  are  inhomo¬ 
geneous  (nonstationary)  Poisson  processes  for  which  the  rate 
function  is  a  random  process.  Such  a  description  is  appropriate 
either  for  the  classical  LV  signal  bursts  with  the  random  turbu¬ 
lent  flow  affecting  the  rate  of  burst  occurrence  or  for  the 
single  photoelectron  pulses  from  the  PMT  with  the  random  classi¬ 
cal  bursts  as  the  rate  function.  Clearly,  when  taken  from  the 
turbulence  to  the  photo-electron  pulses,  a  dual-scatter  LV  signal 
is  a  "triply  stochastic"  filtered  Poisson  process  [11] .  This 
throe  level  nature  of  the  signals  is  illustrated  by  Figure  1. 

Classical  signals  and  burst  counters.  Presently  accepted 
burst-counter  and  frequency-tracker  LV  processors  were  developed 
by  analogy  with  wide-band  frequency  modulation  (FM)  and  radar 
receivers.  For  FM  and  radar  applications  frequency  detection 
(zero-crossing)  circuits  generally  require  about  10  db  signal 
power  to  noise  power  ratio  within  the  bandwidth  of  the  system 
filters.  The  signals  in  such  cases  are  continuous  and  the  noise 

"ijr  ”  ~~  n  -  —  -■  - --  ~  — 

Modified  version  of  noise  model  described  in  [6] . 
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is  stationary  and  Gaussian.  Comments  by  several  speakers  at 
the  1975  Minnesota  LDA  Symposium  indicated  that  they  had 
experimentally  determined  that  their  burst  counter  LV  signal 
processors  failed  significantly  when  the  signal  power  to  noise 
power  ratio  (during  a  burst)  was  less  than  10  db.  This  10  db 
condition  occurs  when  there  are  approximately  10  or  more  photo¬ 
electron  pulses  per  electronic  response  time  x^.  This  response 
time  is  rise  time  or  pulse  width  in  the  case  of  low-pass  fil¬ 
ter.  We  have  also  defined  this  signal  level  as  the  lowest 
value  of  the  "Gaussian"  signal  regime  wherein  the  photomulti¬ 
plier  current  can  be  modeled  as  the  classical  signal  plus  non- 
stat ionary  Gaussian  noise;  (although  neither  Snyder  [9]  or 
Papoulis  [10]  give  any  helpful  rules  as  to  when  this  asymptotic 
approximation  is  valid).  For  lower  signal  power  the  signals 
must  be  treated  as  Poisson. 

Photon  resolved  signals  and  photon  counting.-  A  radically 
different  approach  to  LV  signal  detection  has  been  taken  by 
Pike,  Oliver,  Jakeman,  and  others.  Photon  counting  techniques 
were  developed  for  use  with  low-level  photon  resolved  signals. 
The  summary  results  of  several  years  of  development  of  the 
single-clipping  real-time  photon  correlator  were  described  by 
Oliver  and  Jakeman  in  a  recent  book  [12],  Dr.  Pike  described 
the  application  of  photon-correlation  to  the  processing  of  LV 
signals  at  the  1972  Purdue  conference  [13] .  The  presentation 
was  apparently  not  received  well  by  many  attendees  from  the 
United  States  and  little  has  been  done  in  this  country  with  the 
development  of  photon  counting  techniques  until  recently. 
Increased  interest  was  shown  by  attendees  of  the  1974  Purdue 
Conference. 

One  reason  that  the  single-clipping  correlator  has  been 
slow  to  acceptance  in  this  country  is  that  the  original  theory 
for  its  use  was  based  on  the  assumption  of  many  scatterers  in 
the  probe  volume  with  the  central  limit  theorem  invoked  to 
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render  the  statistics  of  the  scattered  electric  fields  Gaussian. 
Since  this  assumption  is  known  not  to  predominate  in  many  appli¬ 
cations  of  fringe-type  LV  systems  in  air,  none  of  the  first 
theory  for  a  single-clipping  correlator  was  directly  applicable. 
Another  problem  with  the  existing  commercial  correlators  for 
high-speed  air  flow  besides  speed  (minimum  time  resolution  of 
50  nsec)  is  the  lack  of  any  straightforward  way  to  extend  the 
concept  to  the  determination  of  flow  time  statistics  such  as 
correlations  and  power  spectra.  The  system  proposed  in  this 
report  addresses  this  latter  deficiency  as  well  as  eliminating 
the  problems  of  interpretation  of  single-clipping  by  using  full 
multiplication. 

Several  recent  references  provide  additional  valuable 
background  information  on  photon  correlation.  Durrani  and 
Created  [14]  provide  a  derivation  of  the  expected  value  of  a 
photon  correlation  from  single  particle  (Poisson)  signals. 

Birch  et  al.  [15]  have  made  experimental  measurements  in  turbu¬ 
lent  jet  flows  with  skewed  probability  density  functions. 

Abbiss  et  al.  [16]  also  provide  an  analysis  which  shows  that  in 
some  cases  the  Fourier  transform  of  the  correlogram  may  be 
interpreted  as  the  probability  density  of  the  fluid  velocity 
component.’*'  Durrani  and  Greated  [17]  have  investigated  the  use 
of  some  of  the  newer  spectral  estimation  techniques  which  allow 
greater  resolution  from  the  limited  number  of  data  points  in  a 
typical  correlogram.  Finally,  the  reader  should  be  aware  that  a 
new  photon  correlator  instrument  has  been  developed  which  was 
described  by  C.  Fog  [18]  but  the  minimum  time  resolution  inter¬ 
val  is  100  nsec  which  is  rather  slow  for  high  speed  wind  tunnel 
applications.  (It  was  used  for  atmospheric  studies.) 


* 

We  do  not  believe  this  will  be  a  good  approximation  in 
many  practical  cases.  This  is  discussed  in  a  later  section. 
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Burst  rate/amplitude  distributions.-  The  statistical  dis¬ 
tribution  of  ,the  classical  burst  amplitudes  and  the  rate  of 
occurrence  versus  amplitude  are  very  significant  in  the  charac¬ 
terization  of  any  LV  signal  processor.  It  is  generally  accepted, 
for  example,  that  the  optimum  rate  of  occurrence  of  bursts  for 
a  burst  counter  is  less  than  the  inverse  burst  duration  (non- 
overlapping  bursts. )  It  is  also  generally  known  that  the  error 
check  circuits  cause  a  burst  counter  to  emphasize  larger  ampli¬ 
tude  (good  signal  to  noise  ratio)  signals.  On  the  other  hand, 
a  photon  counting  processor  must  emphasize  the  lower  amplitude 
signals  in  a  distribution;  the  higher  amplitude  signals  would 
produce  only  a  single  threshold  crossing  and  otherwise  be 
neglected  by  the  system.  It  is  therefore  not  possible  to  compare 
two  different  types  of  signal  processors  without  knowing  the  sig¬ 
nal  amplitude  distributions  and  the  processor  behavior  as  a 
function  of  burst  amplitude  and  other  factors.  Finally,  in 
order  to  relate  processor  behavior  to  a  specified  particle  size 
distribution,  one  must  first  relate  the  particle  size  distribu¬ 
tion  to  the  burst  amplitude  distributions  and  then  do  all  the 
other  things  already  discussed. 

During  the  initial  phases  of  our  recent  work  for  the  USAF 
Arnold  center  we  addressed  such  questions  as  are  suggested  by 
the  above  statements  both  with  theoretical  models  and  with 
experimental  measurements  of  burst  rate/amplitude  distributions 
for  natural  laboratory  (unfiltered  and  unseeded)  air.  The  reader 
is  referred  to  the  final  report  [1]  for  details.  The  following 
is  one  of  the  concluding  paragraphs  of  that  report: 

The  statistical  distribution  of  the  amplitudes  and 
rates  of  occurrence  of  classical  bursts  has  been  shown 
to  be  central  in  the  problem  of  specifying  or  predicting 
the  data  rates  and  errors  from  any  type  of  LV  signal 
processor.  Differential  and  cumulative  rate/amplitude 
distributions  have  been  formulated  and  analyzed  theo¬ 
retically  and  have  been  measured  experimentally  for  an 
argon  backscatter  LV  system.  The  results  indicate  that, 
for  the  data  obtained,  the  smaller  aerosolc  contribute 
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more  to  the  photon  correlation  accumulator  than  the 
larger  ones.  For  the  data  measured,  there  would 
have  been  available  less  than  300  signals  per  second 
adequate  in  magnitude  to  produce  burst  counter  data 
from  scatterers  larger  than  0.7  ym  in  diameter 
while  there  would  have  been  over  100,000  signals  per 
second  producing  photon  resolved  signals  from  0.2  - 
0.3  pm  diameter  particles. 

Scope 

In  what  follows  we  first  develop  general  Poisson  models 
for  LV  signals  which  include  the  non-st at  ionary  Poisson  occur¬ 
rence  of  photo-electron  pulses  and  the  random  amplitude  effects 
of  both  the  photomultiplier  tube  and  the  particle  scattering 
cross  sections.  Formulas  are  provided  for  conditional  and 
unconditional  moments  including  mean,  variance,  auto-covariance, 
and  higher  moments  (Appendix  A).  These  formulas  are  for  signal 
current,  but  they  become  valid  for  photon  counting  by  a  suit¬ 
able  choice  of  the  PMT  output  filter  impulse  response  function. 

The  next  section  evaluates  the  theoretical  expressions 
for  a  specific  Gaussian  burst  LV  waveform  model.  These  results 
are  used  to  obtain  the  expected  value  of  a  photon  correlation 
estimate.  In  addition  a  Dual  Correlate  and  Subtract  estimator 
which  behaves  as  a  statistical  frequency  discriminator  is  postu¬ 
lated  and  analyzed.  The  following  section  is  devoted  to  statis¬ 
tical  error  analysis  of  the  mean  flow  estimation  technique  using 
the  Dual  Correlate  approach.  The  section  after  that  shows  that 
the  statistical  frequency  discriminator  may  be  applied  to  the 
estimation  of  turbulence  correlations  even  though  the  time  hist¬ 
ory  of  the  velocity  fluctuations  is  not  available  except  as  a 
noisy  randomly  sampled  waveform. 

The  results  of  the  theoretical  considerations  and  the  exper¬ 
ience  we  have  had  previously  with  the  AEDC  [1]  hardware  study 
were  utilized  in  a  system  design  which  is  provided  in  Appendix 
D.  Appendix  C  is  a  derivation  needed  in  the  section  on  variabil¬ 
ity  errors.  Appendix  B  provides  the  theory  and  an  example  pro¬ 
gram  for  correct  Poisson  simulation  of  low-level  LV  signals  for 
evaluation  of  electronic  processor  models. 
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LIST  OF  SYMBOLS 


ENGLISH 


m 

f(  ) 


h 

h(t) 

H(w) 

i 

i(t) 

I 

J 

k 

^(t) 

mie 

M 

aP 

M 

pq 

n(t1(t2) 

nk ' 

N 

Nf 

N. 


b 

P(t) 


bandwidth 

autocovariance  of  i(t) 
electronic  charge 
frequency 

mean  Doppler  frequency 

optical  probe  volume  response  function;  pedastal 
function,  peak  amplitude  normalized  to  unity 

random  single-photoelectron  charge  gain 

Planck's  constant 

impulse  response  of  PMT  and  succeeding  filters 
Fourier  transform  of  h(t) 
counting  index 

photocurrent  at  anode  or  succeeding  filter  output 
optical  intensity 
counting  index 
counting  index 

time  varying  statistical  mean  current  <i(t)> 

nknk-p  "  nknk-q  or  0,ther  function  of  n^ 

photon  correlation  sum  at  delay  PAt 

accumulation  of  (pAT,qAr)  for  Dual  Correlate  Mode 

photoelectron  count  in  interval  (t^,t2) 

photon  count  in  At  interval  about  kAt 

number  of  time  increments  At  (T  *  NAt) 

2 

number  of  fringes  in  the  transmitter  defined  1/e 
probe  volume 

total  number  of  AT  intervals  in  advanced  concept 
operation 

integer  delay  number  in  photon  correlation; 
largest  delay  number  in  dual  correlate  and  subtract 

constant  optical  background  power 

optical  power  incident  on  photocathode 

integer  value  of  delay 
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r,  location  of  the  nearest  approach  of  scatterer  to 

J  center  of  probe  volume 

R  constant  rate  of  signal  bursts 

R(t)  instantaneous  rate  of  signal  bursts 

R^(t)  autocorrelation  of  X(t)  -<X(t  )X(t+T)> 

s(t)  a  generalized  Poisson  shot  noise  signal  (Appendix  A) 

t  time 

t^  instants  of  photoelectron  emission 

T  total  data  collection  time 

inverse  of  mean  Doppler  frequency  l/fm 

u(t)  time  varying  velocity  component 

U(t)  total  velocity  component 

U  mean  velocity  component 

V(t)  vector  velocity 

o 

wQ  1/e  intensity  radius  at  beam  waist 

z(t)  Poisson  Impulse  process 

GREEK 

a  integration  dummy  variable;  also  1/e  half-width 

of  burst 

3  integration  dummy  variable 

6(t )  Dirac  delta  or  unit  impulse  function 

At  simulation  time  resolution  interval 

At  photon  processor  counting  interval 

AT  accumulation  interval  for  second  level  correlation 

e  error 

X(t)  time  varying  mean  rate  of  photoelectron  pulses 

X(t )  -  nP(t)/hv 

Xb  steady  background  photoelectron  rate 

Xj  peak  photoelectron  rate  (pedastal)  from  the  jth 

J  scatterer 

XQ  optical  wavelength 

Xfl  mean  signal  photoelectron  rate  f  <Xj> 

n  product  of  photocathode  quantum  efficiency  and 

dynode  collection  efficiency 
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v  optical  frequency 

p  fractional  turbulence  intensity  o.,/u  *  a  /u 

u  ur  in 

o  rms  deviation  of  the  radian  Doppler  frequencies 

from  scatterers 

2  2 

cA(t)  time  varying  variance  of  i(t):  <(i(t)  -  <i(t)>)  > 

t  delay  variable 

Th  rise  time  or  pulse  width  of  low  pass  h(t) 

Tj  occurrence  time  for  ,jth  signal  burst 

6  beam  intersection  angle;  spherical  declination 

angle  from  direction  of  incident  light 

♦  polar  angle 

U  Poisson  parameter 

u>  Fourier  transform  variable  or  radian  frequency 

random  frequency  of  jth  burst 
um  mean  of  random  variable  utj  (2fffm) 

Special  Notation 

00 

<x>  statistical  expectation  of  x  -  /  xp  (x)dx 

-00  X 

00 

f(t)*g(t)  convolution!  /  f(c)g(t-a)da 

—00 

'  denotes  an  estimate  of  a  statistical  average 
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STATISTICAL  THEORY  OF  DUAL  SCATTER  SIGNALS 
Poisson  Models 

The  Signal  Current.-  The  signal  current  from  a  photomulti¬ 
plier  tube  (PMT)  is  modeled  as  inhomogeneous  filtered  Poisson 
random  process  (see  Appendix  A  and  also  reference  [19] )  given  by 

00 

i(t )  -  l  egjhCt- tA)  (1) 

mm  00 

where  tjL  -  random  time  of  the  ith  photoelectron 
e  -  electronic  charge 
g^  ■  random  charge  gain  of  PMT 
h(t)  ■  impulse  response  of  PMT/filter  system 
The  sy.ivem  response  h(t)  is  obtained  as  a  convolution 
Impulse  response  hp(t),  the  transmission  line  impulse 
h^.(t),  and  the  linear  filter  impulse  response  hf(t) 

h(t)  -  hp(t)*ht(t)*hf(t)  (2)  .. 

where  the  asterisk  denotes  the  convolution  integral: 

f(t)*g(t)  -  /  f (a)g(t  -  a)da  (3) 

.00 

The  superposition  assumes  operation  in  the  linear  range  of  the 
PMT  electron  multiplier.  The  use  of  the  function  hp(t)  assumes 
that  all  single  photo-electron  pulses  have  the  same  shape  except 
for  amplitude.  This  neglects  minor  random  shape  variation. 

The  quantity  which  relates  i(t)  to  the  classical  optical 
power  is  the  statistical  mean  rate  X(t)  of  occurrence  of  the 
randomly  occurring  photoelectron  pulses.  Thus 

X(t)  -  (4) 


of  the  PMT 
response 
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where  n  -  product  of  cathode  quantum  efficiency  and  the 
dynode  collection  efficiency 

hv  ■  Photon  energy 

P(t)  »  Classical  optical  power,  including  background 

light  and  a  constant  component  for  dark  current. 

The  effects  of  dark  current  are  included  by  adding  an 
equivalent  power  Pg.  The  model  could  be  made  more  exact  by  add¬ 
ing  a  separate  dark  pulse  summation  with  a  separate  distribution 
of  amplitudes  which  are  distributed  somewhat  differently  than 
gi;  but  this  distinction  will  not  often  be  required  in  LV  appli¬ 
cations. 

The  previous  material  includes  little  which  restricts  it 
to  LV  signals.  We  now  consider  the  form  of  X(t)  which  is  also 
treated  as  a  filtered  Poisson  process. 

Superposition  of  classical  single  burst  signals.-  Rigorous 
electromagnetic  theory  analysis  of  the  scattered  fields  from 
more  than  one  scatterer  in  the  probe  volume  shows  [20]  mixing 
terms  in  P(t),  the  classical  power  incident  on  the  PMT.  However, 
in  typical  dual-scatter  systems,  the  diffraction  limited  spot 
size  of  the  collecting  lens  is  much  smaller  than  the  probe  vol¬ 
ume;  conservation  of  energy  arguments  show  that  in  such  oases 
the  number  of  scatterers  in  the  probe  volume  may  be  much  greater 
than  unity  with  statistically  negligible  coherent  mixing,  regard¬ 
less  of  the  quality  of  the  collecting  lens  [2] .  This  is  signif¬ 
icant  even  for  LV  systems  which  only  trigger  on  isolated  large 
signal  bursts  because  we  must  also  include  in  the  model  the 
effects  of  smaller  scatterers  which  may  exist  at  higher  number 
density.  We  will  take  the  position  that  at  the  PMT  the  classi¬ 
cal  power  P(t)  is  the  superposition  of  the  backgrouhd  light 
power  and  the  power  from  individual  scatterers  without  coherent 
mixing  cross  terms.  This  will  be  acceptable  so  long  as  the  aver¬ 
age  number  of  scatterers  in  one  diffraction  limited  resolution 
cell  of  the  receiver  is  less  than  unity. 
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A  second  consideration  concerns  the  background  light. 

Even  when  we  neglect  coherent  mixing  of  signals,  there  are  fluc¬ 
tuations  in  the  classical  background  power.  Bertolotti  [21] 
provides  a  review  of  these  effects.  Broadband  background  sources 
can  be  largely  suppressed  by  the  use  of  narrowband  spatial  and 
wave-length  filters,  but  not  always  adequately  enough  for  meas¬ 
urements  from  small  scatterers.  If  the  background  is  modulated 
(for  example  fluorescent  lights)  the  mean  value  signal  is  easily 
removed  by  electronic  filters,  but  the  non-stat ionary  noise  is 
not.  When  the  broadband  background  is  "steady"  there  are  actu¬ 
ally  significant  classical  fluctuations  at  rates  up  to  the  opti¬ 
cal  bandwidth.  Bertolotti  shows  that  when  the  optical  filter 
bandwidth  is  much  greater  than  the  PMT  electronic  bandwidth,  the 
photoelectron  statistics  behave  as  though  the  classical  fluctua¬ 
tions  did  not  exist  (they  are  averaged  out).  Laser  light  scat¬ 
tered  from  windows  is  not  broadband  and  may  exhibit  undesirable 
fluctuations.  This  background  should  be  minimized,  and  its 
effects  studied  further. 

The  random  process  A(t).-  With  cognizance  of  the  preceding 
discussion,  we  model  P(t)  as  the  summation  of  a  constant  back¬ 
ground  Pb  which  Includes  broadband,  laser,  and  dark  current 
sources  and  an  inhomogeneous  filtered  Poisson  signal  process. 

-  p(t>  -  [*b  w]t  <5> 

where  Tj  •  occurrence  time  of  jth  scatterer  reaching  r j , 

Aj  -  random  peak  amplitude  parameter, 

Vj  *  vector  velocity  of  the  jth  scatterer, 

*  location  of  nearest  approach  of  the  scatterer 
J  trajectory  to  the  center  of  the  probe  volume,  and 

f(t,V,r)  *  normalized  optical  system  response  function. 

The  notation  in  equation  (5)  explicitly  shows  that  in  general 

the  shape  of  a  burst  (including  signal  period,  signal  envelope, 
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and  pedestal)  is  a  function  of  the  scatterer  vector  velocity  and 
trajectory  location.  The  response  also  has  random  X^  amplitude 
which  depends  on  both  trajectory  and  particle  scattering  cross 
section.  In  general  the  set  of  instants  Tj  are  independent 
Poisson  random  events  whose  instantaneous  rate,  R(t),  is  statis¬ 
tically  correlated  with  the  velocity  vector. 

Discussion  of  the  model.-  Equation  (5)  is  cast  in  a  gen¬ 
eral  form  which  obscures  certain  details  with  generality.  First 
it  assumes  that  the  velocity  of  a  scatterer  remains  constant 
while  in  the  probe  volume  with  a  value  V ( t j  ,  i* ^  > .  The  extended 
theory  of  filtered  Poisson  processes  is  sufficiently  general  to 
encompass  the  fact  that  the  functional  form  of  the  optical 
response  function  f(t,V,r)  depends  on  two  vector  random  varia¬ 
bles.  *  However,  Snyder  [9]  assumes  that  the  vector  random 
parameters  are  independent.  We  are  not  certain  at  the  present 
time  what  the  statistical  dependence  of  the  rate  function  R(t) 
on  the  velocity  ?(t)  may  imply,  but  no  serious  consequences  will 
result  with  low  turbulence  flow, 

Conditional  Signal  Statistics  of  the  Photocurrent 

At  times  the  models  for  the  systems  analysis  problem  may 
be  simplified  until  analytical  methods  are  applicable.  In  these 
cases  the  use  of  conditional  statistics  will  usually  simplify 
the  analysis.  Papoulls  [10]  discusses  the  use  of  conditional 
statistics  at  length.  We  utilize  thiB  technique  at  length  in  a 
later  section.  Basically  for  a  multilevel  random  process  the 
technique  consists  of  assuming  the  higher  level  random  processes 
are  known  and  deterministic,  evaluating  conditional  expectations 
assuming  the  higher  level  processes,  then  evaluating  the  expec¬ 
tation  of  the  result  with  respect  to  the  higher  level  processes. 
First  we  will  consider  statistics  of  i(t)  assuming  the  classical 
optical  signal  X(t)  is  known. 

♦Elementary  shot  noise  theory  is  restricted  to  an  impulse 
response  function  which  is  constant  in  shape. 
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Instantaneous  mean,  variance,  autocovariance.-  The  result 
In  Appendix  A  may  be  applied  to  determine  the  instantaneous  mean, 
variance  and  auto-covariance  of  the  signal.  These  are  given  in 
terms  of  the  function  X(t).  The  results  are  as  follow: 

mt(t)  -  <i(t)>  -  e<gi>X(t)*h(t)  (6) 

c12(t)  -  <(i(t)  -  <i(t)>)2>  -  e2<g2  >X(t)*h2(t)  (7) 

cii(trt2)  -  <i(t1)i(t2)>  -  <i(t1)><i(t2)>  (8) 

00 

■  e2<gi2  >  /  X(a)h(t1  -  a)h(t2  -  a)da 

where  <  >  denotes  statistical  expectation  and  where  the  asterisk 
again  denotes  the  convolution  integral.  These  results  include 
the  specification  that  hp(t),  the  impulse  response  of  the  PMT 
anode,  have  unit  weight,  i.e., 

!  h  (t)dt  -  1  (9) 

-00  “ 

in  order  to  maintain  conservation  of  charge.  The  functions 
ht(t)  and  hf(t)  may  include  amplification  or  loss  factors  and 
need  not  have  unity  weight. 

Conditional  noise  and  SNR.-  The  concept  of  signal-to-noise 
ratio  arose  in  communications  theory  when  the  "noise"  was  an 
additive  stationary  Gaussian  random  process  totally  character¬ 
ized  by  a  mean,  mean-square  deviation  (variance),  and  a  power 
spectral  density.  The  ratio  of  the  peak  or  average  signal  to 
the  rms  noise  was  a  useful  measure.  The  preceding  equations 
show  the  mean-square  deviation  (variance)  to  be  an  instantaneous 
time  function  which  is  related  to  the  classical  signal.  Observa¬ 
tion  of  real  LV  signals  on  an  oscilloscope  display  or  computer 


IB 


simulations  such  as  that  shown  in  Figure  2  show  that  the  concept 
of  slgnal-to-noise-ratio  is  not  an  adequate  figure  of  merit 
without  careful  specification. 

For  an  example  SNR  definition,  we  consider  a  low-pass  PMT 
impulse  response  as  a  rectangular  function: 


h(t)  -  -r-  Rect  (t/T.  )  (10) 

Th  n 


where 


Rect  (t)  B  1,  1 1 1  _<  0 . 5  (11) 

-  0,  | t |  >  0.5 


If  we  now  also  assume  that  tu  <<  TL  where  T  is  the  signal 

h  m  m 


Figure  2.  Computer  Simulation  of  LV  Signal 

using  Algorithms  similar  to  Appendix 
B.  (By  J.  F.  Meyers  NASA  Langley.) 
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period  of  interest,  then  we  could  obtain  an  instantaneous  SNR 
from  equations  (6)  and  (7)  as 


SNR(t)  »  m2  (t)/o2  (t) 


<g1>2X2(t)/<g1^>  X(t)/x, 


(12) 


<8±> 

*  — : X(t)T, 

<er>  1 


2  2  * 
where  the  quantity  <g^>  /<gi  >  is  typically  between  0.5’  and 

1.0  with  magnitude  depending  on  the  relative  variance  of  the 

PMT  single  photoelectron  pulse  gain.  For  an  ideal  tube  the 

quantity  X(t)xh  would  be  the  Instantaneous  SNR.  This  is  not 

useful  since  it  is  a  time  function  instead  of  a  number. 

As  an  alternative,  we  may  take  the  local  time  average  of 
the  SNR  given  by  (12)  over  a  single  cycle  near  the  peak  of  the 
pedastal.  This  would  give,  for  an  ideal  PMT, 


SNRAVpeak  XjTh  (13) 

where  Xj  is  the  peak  value  of  the  pedastal  of  the  jth  signal 
burst,  if  we  assume  sparse  non-overlapping  bursts.  We  observe 
that  equations  (6)  through  (8)  are  valid  when  h(t)  is  a  bandpass 
function,  but  (13)  is  then  meaningless  unless  we  redefine 
for  a  bandpass  h(t).  Also  we  note  that  this  definition  would 
be  necessary  for  meaningful  use  with  a  burst-counter  processor, 
since  it  is  the  bandpass  filtered  AC  signal  to  wlde-band  noise 
power  that  is  significant  in  that  case. 

Signal  Regimes.-  The  idealized  quantity  SNRav  peak  8lven 
in  equation  (13)  as  Xjth  is  at  least  a  useful  quantity  in  defin¬ 
ing  a  classification  system  of  signal  regimes  for  a  low-pass 


* 


Typical  rms  values  are  0.707  or  greater  ([22], page  86). 
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filtered  signal.  The  following  definitions  of  a  photon-resolved, 
a  photon- limited,  and  a  Gaussian  signal  regime  have  been  some¬ 
what  arbitrarily*  Identified: 

The  signal  is  photon  resolved  if  Xt^  <<  1.  In  this  case 
the  probability  of  two  or  more  photoelectron  events  occurring 
within  the  response  time  Th  is  small.  Its  appearance  is  that 
of  individual  pulses  which  vary  in  height  due  to  the  randomness 
of  g^.  Photon  counting  methods  are  appropriate.  The  condi¬ 
tional  mean  value  of  i(t)  is  still  proportional  to  X(t),  but 
there  is  no  visible  resemblance  to  X(t).  This  condition  is 
illustrated  by  the  extreme  right  hand  portion  of  Figure  2. 

For  Xth  >>  1,  the  signal  i(t)  is  asymptotically  a  nonsta¬ 
tionary  Gaussian  Process.  In  this  case  the  first  and  higher 
order  probability  density  functions  for  i(t)  at  any  set  of 
instants  ('tiito*  ♦  •  • )  may  be  determined  immediately  by  plugging 
the  mean,  variance  and  auto-covariance  from  the  preceding  equa¬ 
tions  into  well  known  Gaussian  formulas.  Under  these  same  con¬ 
ditions  the  signal  display  appears  to  the  eye  as  a  classical 
signal  m^t)  plus  Gaussian  noise.  This  condition  is  approached 
by  peak  of  the  trace  in  Figure  2.  The  major  difference  between 
this  case  and  that  of  classical  communications  theory  problems 
is  that  the  a  value  for  the  noise  is  signal  (time)  and  system 
dependent.  Usually,  signals  in  the  Gaussian  regime  are  suitable 
for  processing  by  classical  methods  (burst  counter  and/or 
tracker) . 

The  photon  limited  regime  is  that  for  which  Xrh  is  within 
an  order  of  magnitude  of  unity.  No  mathematical  simplifications 
are  possible.  Visually  the  signal  appears  as  shown  in  all  but 
the  lowest  portions  of  Figure  2.  The  upper  limits  of  photon 


Gee  Papoulis  [10]  page  571,  No  specific  limit  on  the 
magnitude  of  Xt.  is  given.  We  define  the  photon  limited  regime 
as  0.1  <  Xth  <10. 
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counting  techniques  and  the  lower  limits  of  conventional  tech¬ 
niques  both  fall  in  this  range. 

As  we  have  illustrated,  the  signal  classification  may 
apply  to  different  portions  of  the  same  waveform.  We  may  also 
use  the  classification  with  respect  to  the  peak  pedastal  value 
XjTfa  to  classify  signal  bursts,  Under  this  type  of  classifi¬ 
cation,  Figure  2  illustrates  a  photon-limited  burst  whose  peak 
average  SNR  is  less  than  10.  Additional  bandpass  filtering 
would  increase  and  place  the  central  portion  of  the  burst  in 
Figure  2  in  the  Gaussian  regime.  This  would  not  be  possible 
with  a  significantly  weaker  signal. 

Unconditional  Statistics 

Long  time  mean,  variance,  autocovariance.-  Equations  (6)  - 
(8)  include  the  assumed  deterministic  classical  signal  X(t)  which 
is  proportional  to  instantaneous  optical  power.  When  we  wish 
later  to  evaluate  the  long-time  average  result  which  accumulates 
during  a  photon  counting  experiment,  it  will  be  necessary  to 
treat  X ( t )  as  an  ergodic  random  process  with  long-time  average 
equal  to  the  unconditioned  statistical  mean: 

<X( t )>  -  X  (14) 

We  also  make  use  of  the  autocorrelation  of  X(t): 

<X(t)  X(t  +  t)>  *  Rx(t)  (15) 

Now  from  equation  (6)  taking  the  expectation  with  respect  to 
X ( t )  gives  the  average  current  as 

00 

<i>  ■  e<gi><X>  /  h(t)dt  (16) 
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where  the  integral  is  unity  unless  h(t)  includes  preamplifi¬ 
cation  or  attenuation  external  to  the  PMT.  In  order  to  deter- 

2 

mine  the  long  time  variance  a.  we  do  not  get  the  correct 

1  2 

answer  by  taking  the  expectation  of  o,  (t)  given  by  equation 

1  2 

(7).  Rather,  one  determines  the  conditional  value  of  <i  (t)> 
by  adding  the  square  of  (6)  to  (7).  The  expectation  with 
respect  to  X(t)  follows;  finally,  the  square  of  equation  (16) 

p 

is  subtracted  from  the  unconditional  expectation  of  i  (t).  When 
all  these  steps  are  completed,  and  similar  ones  for  the  uncondi¬ 
tional  autocovariance,  we  obtain: 

non00  a  00  o 

V  -  eZ[  <gi>2  /  Rx(a)fh(ct)da  +  <8*><\>  /  hz(a)da  (17) 

—  00  —00 

00 

-  <gt>2<x>2  (/  h(cOda)2] 

4  1  —00 

Cii(T)  -  e2[<gi>2Rx(T)*fh(T)  +  <gi2><X>fh(T)  (18) 

00 

-  <g  1>2<X>2(/  h(a)dct)2] 


where 


00 

fh(T)  "  /  h(a)h(a + r )da  (19) 

—00 

The  second  term  in  the  expression  for  ^^(t)  vanishes  for  t 
greater  than  the  impulse  response  time  for  the  PMT  and  filter 
combination;  the  last  term  is  the  square  of  the  mean;  the  first 
term  is  the  correlation  of  X(t)  smoothed  by  Ihe  correlation  of 
h(t )  with  itself. 

Ideal  photon  correlation.-  An  idealized  photon  correlator 
counts  all  photoelectron  emission  events  during  successive  uni¬ 
formly  spaced  clock  periods  of  duration  At.  The  number  sequence 
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{n^}  which  result  _s  algebraically  manipulated  to  yield  the 
summation  of  terms  nknk+p.  In  evaluating  the  expected  value 
of  the  result  of  accumulating  such  a  sum  we  encounter  the  need 
to  evaluate  the  quantities  <nk>  and  <nk  +  nk+p>.  These  expec¬ 
tations  may  be  evaluated  using  equations  (16)  -  (19)  by  assuming 

h(t )  -  Rect (t/Ax )  (20) 

■  1/e 

where  Rect(t)  was  defined  in  equation  (11)  and 

/ h(t)dt  ■  At  (21) 

With  these  assumptions,  i(t)  is  equal  to  the  number  of  photo¬ 
electron  events  in  the  interval  (t- At/2,  t  +  Ax/2)  and  the 
formulas  reduce  to 


<i>  ■  <n^>  ■  <A>At 


00 

var  nk  -  a  2  -  <n^  -<nfc>2  *  At  /  Rx(ct)  A(^)da 


+  <A>At  -  <A>2At2 


(22) 

(23) 


<nknk+P>  “  cii‘p4T)  +  <1>2  <24> 

-  ATftx(T)»A(T/&T)|T.piTp)l0 

where  the  correlation  integral  of  equation  (19)  produces  a  tri¬ 
angular  function,  i.e. 

/  Reot  (^)  Rect(^I)da  -  AtA(t/At)  (25) 
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where 


A(t/At  )  -  1  -  |  t  |/  At  ,  |  t  |  <  At  (26) 

-  0  ,  |  t  |  >  At 

We  observe  that  the  generally  accepted  result  that  the  photon 
correlation  Is  shaped  like  the  correlation  of  the  classical 
signal  Is  true  subject  to  the  jump  discontinuity  at  zero  delay 
and  subject  to  the  triangular  weighting  function  which  behaves 
as  a  low-pass  filter  with  respect  to  the  details  of  R^(t). 

When  At  is  much  smaller  than  a  characteristic  signal  per¬ 
iod,  then  equations  (23)  and  (24)  simplify: 

<nknk+p>  “  At2  rx<P*t>>  P  t  0  (27) 

«  At2  R^(0)  +  <X>At,  p*0 

var  n^  *  ■  <X>At  +  At2  (<X2>  -  <X>2)  (28) 


i.e.  , 


o 

var  nfe  -  <n^>  +  At  var  X  (20) 

This  last  result,  which  we  obtain  as  a  special  case,  has  also 
been  given  by  Bertolotti  [21] .  It  provides  a  way  to  measure  the 
variance  of  the  classical  signal  even  with  photon  resolved  sig¬ 
nals  if  a  long  sequence  of  nk  values  are  available. 

PHOTON  COUNTING  PROCESSORS  FOR  MEAN 
FLOW  AND  TURBULENCE  INTENSITY 

In  this  section  we  provide  an  idealized  theoretical  basis 
for  the  use  of  photon  correlation  and  a  new  type  of  photon 
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counting  technique,  "dual-correlate-and-subtract , "  and  show  how 
they  may  be  used  for  mean-flow  and  turbulence  Intensity  measure¬ 
ments.  Extensions  to  higher  order  statistics  are  discussed  in  a 
later  section.  Error  analysis  of  the  processors  described  in 
this  section  is  presented  in  the  next  section. 

Specific  Signal  Model 

The  form  of  equation  (5)  is  quite  general.  In  most  prac¬ 
tical  systems,  the  function  f(t,?j,Fj)  which  describes  the  opti¬ 
cal  response  with  respect  to  particle  position  and  velocity  is 
complicated  when  the  effects  of  limiting  pinhole  apertures  and 
variable  duration  due  to  high  turbulence  are  included.  For  the 
present  we  assume  a  simplified  low  turbulence  model  which  assumes 
a  burst  with  perfect  contrast  and  constant  shape: 

X(t>  -  Xb  +  l  Xjf(t  -  Tj)[l  +  COS  U)j(t-Tj)]  (30) 

where  Xj  -  random  burst  pedastal  amplitude 

Xfe  ■  constant  background  rate, 

f(t)  -  low-pass  pulse  waveform  with  peak  equal  unity, 

tj  -  occurrence  time  for  jth  burst, 

w.  «  radian  frequency  proportional  to  one  velocity 
J  component  of  the  jth  particle 

We  may  write  equation  (30)  as  the  sum  of  a  constant,  X^,  a  low- 
pass  process,  X^^(t),  (the  pedastals),  and  a  bandpass  process 
XbpOO.  Then 

X(t)  -  Xb  +  X£p(t)  +  Xbp(t)  (31) 

We  assume  that  there  are  several  fringes  in  the  probe  volume  so 
that  the  spectra  of  X^p(t)  and  Xbp(t)  are  non-overlapping.  Thus 
Xbp(t)  is  a  zero-mean  process,  and  XAp(t)  and  Xbp(t)  are  uncor¬ 
related.  We  obtain 
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(32) 


BX(T)  -  <*(t)>2  +  CUp(T)  +  CXbp(T) 

We  now  use  the  low  turbulence  assumption  and  further  assume 
that  the  scatterers  are  uniformly  dispersed  In  space  so  that, 
R(t),  the  rate  of  burst  arrivals  Is  a  constant  R.  The  results 
of  Appendix  A  can  be  applied  to  derive  expressions  for  the 
three  terms  In  equation  (32),  the  result  Is 

00 

<X(t)>  -  Xb  +  R <X . >  /  f (t)dt  (33) 

w  m  00 

CUp(T)  '  B<*J>  /  *<*)  1  (34) 

00 

CXbp^  "  R<^j>  i  <oos  WT>  /  f(t)f(t  +  T)dt  (35) 

where  the  expectation  of  the  <cos  wt>  term  is  with  respect  to 
the  random  variable  u>j.  The  derivation  requires  that  we  expand 
the  product  of  cosines  with  the  sum  and  difference  formula  and 
approximate  the  integral  of  the  product  of  a  low-pass  term  and 
a  bandpass  term  as  zero. 

We  now  assume  that  the  turbulence  is  Gaussian  with  mean 
radian  frequency  u>m  and  rms  deviation  a  .  Then  by  direct  appli¬ 
cation  of  the  definition  of  expectation  we  obtain41 

os  -(id— w  )2/2a2 

<cos  wt>  »  / - e  m  w  cos  a)T  dw  (36) 

*  /IS?  cr 

0) 


■  e  cos 

m 


*In  general  <cos  wt>  is  simply  related  to  the  statistical 
characteristic  function  for  the  random  variable  u.  since  this 
function  is  defined  as  a  Fourier  transform  of  theJprobabillty 
density  function  [10] . 
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A  simplified  expression  for  the  autocorrelation  of  the  classical 

signal  is  thus 


Rx(t)  -  <Xb  +  R<Xj>  /  f(t)dt)2  (37) 

n  -j  -<J2  T2/2  « 

+  R<X“>(1  +  I  e  w  cos  wmT )/  f(t)f(t+x)dt 

This  result  shows  that  regardless  of  the  shape  of  the  envelope 
function  f(t),  the  autocorrelation  function  has  a  consinusoidal 
variation  at  the  mean  signal  frequency  and  a  Gaussian  envelope 
decay  factor  due  to  turbulence  Intensity. 

If  we  now  assume  that  the  classical  bursts  are  Gaussian 
shaped  (TEMo0  beams  without  aperture  effects),  then  we  have 

2  2 

f<t)  -  e-t  I*  (38) 

where  a  is  the  1/e  half  width  of  the  envelope  and  obtain 

00 

/  f(t)dt  ■  /H  o  (30) 

.00 

/  f(t)f(t + x)dt  -  ae“t2/2“2 

The  final  simplified  expressions  for  the  first  two  moments  of 
X(t )  are 


<X>  »  <X(t)>  ■  Xb  +  R*Xj>  /Tr  a 

2  2 


Rx(t)  -  <X>2 +  R<X2>^Ja(l  +  |  e  ^cos  u>mT)e“t2/2a2 


where  low  percent  turbulence  has  been  assumed  and 


(40) 

(41) 


25 


Ab  ■  background  count  rate, 

R  -  scatterer  arrival  rate, 

Xj  ■  random  pedaatal  height  from  jth  scatterer, 
a  -  1/e  half  burst  duration, 

■  rms  deviation  of  radian  frequency  due  to  turbulence, 
wm  *  mean  radian  frequency  «  1),  and 

t  e  delay  variable  of  autocorrelation. 

The  shapes  of  typical  correlation  functions  for  zero  tur¬ 
bulence  and  10%  turbulence  are  illustrated  in  Figure  3. 


Idealized  Photon  Correlation  of  LV  Signals 

The  number  n^  is  the  number  of  photo-electron  emission 
events  in  the  interval  which  extends  ±At/2  from  the 
instant  kAt.  An  idealized  photon  correlator  produoes  and 
sums  delayed  products  from  the  uniformly  spaced  sequence  (n^h 
Ve  assume  the  total  number  of  products  accumulated  is  N.  The 

A 

accumulator  produces  a  sum  Mp  at  the  delay  value  pAt  given  by 

fip  ■  X  Vk+p  <42) 

The  ideal  photon  correlator  would  simultaneously  accumulate  Nnk 
defined  by 

A  N 

Nnk  "  j*  nk  (43) 

The  unconditional  expected  value  of  these  sums  is  obtained  from 
equations  (22)  and  (27)  after  interchanging  expectation  and  sum 
mat ion  as 


<M  >  •  N<n,  n,,  > 
p  k  k+p 


-  NAT2Rx(pAT),  p  f  0 

-  N[<\>At  +  At2Rx(0)] ,  p 


0 


(44) 
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Turbulence  ■'  10% 


(45) 


<Nnk>  =  N<>>At 

By  combining  equations  (40)  -  (45),  we  may  obtain  an  esti¬ 
mator  for  the  autocorrelation  of  >.(t),  i.e., 


ft.(pAx)  -  - i- 

A  AtzN 


1 

At*N  k 


N 

I  nkn 

-1  * 


k+p' 


P  f  0 


(40) 


&x(PAt) 


1 

-75 

AtzN 


(fllp  -  Nnk),  p  -  0 


1 


At*N 


I  MV"1)*  p  “  0 

k-1  K  K 


This  estimator  includes  the  zero  delay  value,  which  is  usually 
emitted,  by  making  use  of  the  separate  mean  count  computation. 
The  same  mean  count  estimate  nk  may  be  used  to  estimate  the  long 
delay  level 

<nk/At>2  -  <X(t)>2  -  Rx(«)  (47) 


Interpretation  of  an  autocorrelation  estimate  computed 
according  to  equation  (46)  involves,  first,  the  use  of  an  ana¬ 
lytical  model  such  as  equation  (41);  second,  a  parameter  extrac¬ 
tion  procedure,  such  as  a  mean  square  error  minimizing  curve  fit 
algorithm;  third,  a  correction  for  any  statistical  bias  errors; 
and  fourth,  a  variability  error  criterion  which  assures  that 
sufficient  data  is  accumulated.  We  have  provided  a  procedure 
and  an  example  model  for  the  first  step.  The  bias  and  varia¬ 
bility  errors  are  discussed  in  a  later  section.  We  have  not 
considered  the  optimization  of  the  second  step  although  one 
method  is  discussed  below.  Some  literature  [15,10,17] is  begin¬ 
ning  to  appear,  but  further  effort  is  needed. 
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At  the  Denmark  LDA  conference  in  the  summer  of  1975, 

Abies  et  al.  [16]  described  a  "new"  interpretation  of  a  photon 
correlellogram  as  the  Fourier  transform  of  the  velocity  proba¬ 
bility  density  function.  This  "breakthrough"  allows  the  meas¬ 
urement  of  the  probability  density  of  the  velocity  field  by 
digital  Fourier  transform  of  the  correlation  results  under  the 
condition  that  the  pedestal  and  fringe  envelope  correlations 
are  nearly  constant  over  the  range  of  delays  for  which  the  sin¬ 
usoidal  correlation  is  appreciable.  This  result  is  equivalent * 
to  the  interpretation  of  spectrum  analyser  displays  (connected 
directly  to  the  photo-detector)  years  ago  as  the  probability 
density  of  the  velocity  field.  The  restriction  is  equivalent 
to  requiring  many  fringes  in  the  probe  volume  so  that  the  transit¬ 
time  spectral  broadening  is  small.  Our  equation  (35)  shows  the 
relationship  between  the  probability  density  function  of  the 
velocity  (frequency)  samples  and  the  correlation.  The  expres¬ 
sion  <oos  wt>  is  the  cosine  Fourier  transform  of  the  probability 
density  for  When  there  are  many  fringes  in  the  probe  volume 

the  envelope  correlation  term  is  a  broad  pulse  which  is  essen¬ 
tially  constant  over  the  extent  of  <cos  wt>.  In  that  case,  the 
transform  of  the  correlation  is  actually  inverting  the  statis¬ 
tical  characteristic  function  tp  produce  a  scaled  probability 
density. 

The  unfortunate  truth  remains  that  in  most  practical  LV 
problems,  the  probe  volume  and  spatial  frequency  of  the  fringes 
must  be  small  with  the  result  that  we  rarely  have  the  luxury 
of  having  many  fringes  in  the  probe  volume.  For  the  realistic 
case,  then,  we  find  the  spectrum  to  be  the  convolution  of  the 
velocity  probability  density  function  and  the  transform  of  the 
envelope  correlation,  and  care  is  required  in  interpretation. 

- — — - 

At  least  in  theory:  the  photon  correlator  is  much  more 
efficient  at  low  signal  levels  than  any  swept  frequency  spectrum 
analyser  could  be. 
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A  Photon  Counting  Frequency  Discriminator 

This  section  describes  the  basis  for  a  statistical  mean 
frequency  discriminator  for  photon  resolved  LV  signals.  We 
begin  with  a  historical  description  of  the  motivating  logic 
which  may  offer  insight  to  others  for  more  advanced  development. 

Autocorrelation  Frequency  Discrimination.-  It  is  well 
known  that  the  autocorrelation  function  of  a  zero-mean  narrow- 
band  random  process  is  cosinusoidal  with  frequency  equal  to 
the  mean  frequency  of  the  narrow-band  process  and  envelope 
which  is  the  autocorrelation  of  the  amplitude  and  phase  envelope 
of  the  process.  For  example,  if  the  two-sided  power  Sv(f)  of 
a  random  process  x(t)  is 

Sx(f)  -  A(f-fm)  +  A(f  +  fm)  (48) 

Where  A(f)  is  an  even  pulse-shaped  function  and  fm  is  a  center 
frequency  much  greater  than  the  frequency  width  of  A(f),  then 

Rx(t)  -  2Ra(t)  cos  27rfmT  (40) 

Where  A(f)  is  the  Fourier  transform  of  Ra(x).  Now  we  note  that 

Cl 

for  small  perturbations  of  the  argument  2irxfm  about  the  point 

t  -  3/4  fm  (BO) 

the  value  of  Rx(x)  varies  in  proportion  to  the  perturbation  of 

either  x  or  fm<  Thus,  if  x  is  selected  to  satisfy  equation 

(80),  then  the  statistical  autocorrelation  function  would  serve 
Instantaneously*  as  a  frequency  discriminator  for  small  deviations 

4  1 

This  is  significant  in  later  sections  even  though  the 
statistical  autocorrelation  function  cannot  be  experimentally 
observed. 
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of  fm.  If  the  random  process  were  ergodic,  then  the  same 
results  would  apply  to  time-average  autocorrelations  when  the 
averaging  time  was  long  compared  with  the  coherence  time  of  the 
process.  These  principles  may  be  extended  to  a  random  process 
which  is  the  sum  of  a  narrow-band  process  and  a  low-pass  pro¬ 
cess  by  first  filtering  the  process  with  a  high-pass  filter  to 
remove  the  non-zero  mean  low-pass  process. 

Application  to  photon  resolved  signals.-  In  the  case  of 
photon-resolved  LV  signals,  the  classical  signal  random  process 
is  not  necessarily  recoverable  from  the  sparse  single  photon 
events,  but  we  have  already  shown  that  the  autocorrelation 
function  may  be  approximated  by  the  expected  value  of  a  photon 
correlation  operation  on  the  photon  events.  The  original  con¬ 
cept  was  thus  to  devise  a  high-pass  digital  filter  which  would 
be  applied  to  the  count  sequence  {n^}  to  remove  the  effects  of 
background  light  and  low-pass  pedastal  from  the  statistical 
average  <nk>  while  leaving  the  bandpass  information  signal 
information.  This  filter  would  be  followed  by  a  digital  corre¬ 
lation  (delayed  product  summation)  at  the  delay  value  given  by 
equation  (SO).  Since  the  digital  electronics  had  to  be  very 
fast,  only  simple  digital  filters  could  be  considered.  The 
simplest  one  that  we  thought  of  was  to  delay  the  sequence  {n^} 
by  one-half  period  of  the  signal  and  subtract.  This  operation 
cancels  the  low-frequency  portions  of  the  expected  value  of 
nk  and  adds  the  sinusoidal  portions  in  phase.  Such  a  delay  and 
subtract  filter  would  produce 

mk^  "  nk  ”  nk-q  (81) 

Where  the  counting  interval  At  must  be  adjusted  to  satisfy 

qAt  -  l/(2fm)  (52) 
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with  q  an  even  Integer.  After  this  was  done  one  would  accumu¬ 
late  lag  products  of  m^,  with  delay  pAi  equal  to  three  quarters 
of  a  signal  period  p  -  3q/2.  This  approach  leads  one  to  form 
the  summation  M'pq  given  by 


ft 


pq 


X  “  nk-pnk-q  +  2nknk-p  -  Vk-p-q 


(53) 


The  previous  results  for  photon  correlation  are  applicable  for 
the  determination  of  the  expected  value  of  the  quantity  ft'pq. 

The  details  are  omitted  here  and  included  for  the  simpler  dual¬ 
correlate  approach  described  below,  but  Figure  4  is  a  simplified 
discriminator  characteristic  for  the  expected  value  of  the 
result.  It  is  provided  for  comparison  with  Figure  5  which  is 
described  in  detail  below. 


After  all  of  the  above  reasoning,  it  occurred  to  us  that 
another  approach  (we  thought)  consisted  of  performing  the  delay 
and  subtract  filtering  operation  on  the  correlation  estimate 
after  it  was  made  rather  than  on  the  high-speed  signal  sequence 
{n^}.  Reference  to  Figure  3  shows  that  a  one-half-period  shift 
and  subtract  of  the  typical  LV  correlation  will  approximately 
cancel  the  low-pass  portion  of  the  autocorrelation.  The  results 
of  this  approach  provide  an  approximate  discriminator  response 
as  shown  in  Figure  8  with  less  high-speed  data  arithmetic 
required.  This  "dual  correlate  and  subtract"  approach  also 
leads  to  larger  frequency  range  for  the  discriminator  function, 
and  thus  has  been  chosen  for  development.  Hindsight  has  now 
shown  us  that  by  algebraic  rearrangement  this  approach  is  most 
economically  implemented  with  high-speed  delay  and  subtract 
filter  prior  to  multiplication.  This  technique  may  be  shown 
to  be  just  another  implementation  of  the  original  filter  and 
correlate  concept. 
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Dual  correlate  and  subtract.-  We  will  let  the  quantity 

be  defined  as  follows 
P4 


pq 


N 

“  I 

k-1 


m. 


CL^-p)  ■  ( I0*”*-* 


(54) 

•1 

« 


where 


“k  ‘  Vnk-p  -  “k-q1  <B6> 

It  is  straight  forward  to  show  that,  except  for  a  few  end 
terma'Aq  thematically  identical  to  the  quantity  Mp  -  Mq 
where  and  are  defined  by  equation  (42).  It  is  for  this 
reason  that  we  will  lable  the  approach  as  the  "Dual  Correlate 
and  Subtract"  technique  even  though  the  delay  values  are  nega¬ 
tive  instead  of  positive.*  We  now  demonstrate  that  the  expected 
value  of  Mpq  behaves  as  a  frequency  discriminator  as  illustrated 
in  Figure  5  under  conditions  which  we  will  identify!  The 
adjustment  of  the  system  clock  period  At  leads  to  a  null  in  the 
expected  accumulator  value.  Measurement  of  At  provides  a  direct 
measure  of  the  mean  signal  frequency  as  we  shall  now  show  with 
our  simplified  signal  models. 

From  equation  (44)  we  obtain  the  expected  value  of  the 
estimator  as 

<flpq>  -  NAt2[Rx(PAt)  -  Rx(qAT)J  (58) 

The  complete  expression  is  obtained  for  our  simplified  signal 


The  negative  delay  implementation  was  more  suitable  for 
the  hardware  design. 
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model  by  using  equation  (41)  in  its  entirety.  Here  we  assume 
that  pAt  and  qAc  are  both  small  compared  with  a  so  that  the 
pedastal  terms  cancel  as  well  as  the  steady  term.  This  leaves 

<fl  >  -  N4tZR<X*>  a  .  (57) 

-  3  <V°S>(pat>2  - 

•  e  a  cos  w,j?At  -  e  a  cos  w  aAx  . 

m  m 


Now  we  require  that 


p  ■  3q  (58) 

and  observe  the  behavioi  of  equation  (57)  near  the  values  of 
At  specified  by  letting  qAt  be  one  quarter  of  the  signal  per¬ 
iod  where  both  cosine  terms  will  vanish. 

qAx  «  27T/4uim  -  Tm/4  (59) 

The  shape  of  the  term  in  braces  is  plotted  in  Figure  5 
under  the  assumption  of  many  fringes  in  the  probe  volume 
(a  large)  and  low  turbulence  (o^  small).  Thus  the  quantity 
plotted  is  simply 


[cos  3qAt  <t>  -  cos  qAt  w  1  (60) 

m  m 

Figure  5  illustrates  the  expected  behavior  of  the  accumulator 
sum  for  changes  in  the  mean  signal  frequency  t*)m.  If  the  sys¬ 
tem  clock  frequency  is  changed  to  change  At ,  then  the  response 

2 

is  the  product  of  At  and  the  curve  shown  in  the  figure.  The 
shape  of  the  curve  is  affected  but  the  zero  crossing  locations 
are  not . 
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Selection  of  delay  constants.-  The  theory  does  not 
uniquely  specify  the  relationship  of  At  (the  above  system  clock 
period)  to  Tm  the  signal  period  because  q  is  not  specified. 

For  a  given  value  of  signal  period,  T  ,  the  largest  possible 

j.  m 

value  of  At  for  an  acceptable  null  is  when  q  -  1  and 
At  ■  Tm/4.  This  produces  the  least  variability  error  (as  we 
will  show)  and  the  most  bias  error  due  to  time  smear  (the  tri¬ 
angular  weighting  function  in  equation  (24).  The  value  q  «  1 
also  allows  the  highest  signal  frequencies  to  be  measured  for 
a.,  given  maximum  system  clock  frequency.  The  bias  errors  can 
be  reduced  at  the  expense  of  increased  variability  error  and 
reduced  maximum  signal  frequency  by  using  q  ■  2,  At  -  Tm/8. 

In  brdcr  to  facilitate  experimental  research,  our  design  for  a 
research  Instrument  allows  selection  of  p  and  q  over  a  wide 
range. 


STATISTICAL  ERROR  ANALYSIS 
Introduction 

The  two  principle  types  of  error  which  arise  in  statis¬ 
tical  measurements  are  bias  error  and  variability  error.  Bias 
error  is  a  term  which  refers  to  the  difference  between  the 
statistical  expectation  of  the  measurement  system  output  and 
the  desired  average  value  being  measured.  The  variability 
error  is  the  rms  value  of  the  random  deviation  of  a  specific 
experimental  result  from  the  statistically  expected  value.  For 
ergodic  random  processes,  the  variability  error  converges  to 
zero  in  the  limit  of  infinitely  long  data  collection  time;  but 
it  converges  to  an  acceptable  level  (which  must  be  defined  in 
the  measurement  objectives)  within  a  finite  measurement  time. 
The  bias  error  cannot  be  removed  by  further  averaging  but  it 
can  often  be  removed  by  analytical  compensation  or  by  experi¬ 
mental  calibration  when  it  is  small  compared  with  the  desired 


36 


quantity.  In  general,  analysis  of  both  types  of  errors  is 
required  in  any  statistical  measurement. 

There  are  three  different  methods  for  evaluation  of  the 
statistical  errors  of  a  measurement  system;  these  are:  analysis, 
simulation,  and  experiment.  Appendix  A  provides  theory  for  the 
higher  order  moments  of  inhomogeneous  Poisson  processes  which 
we  have  expanded  for  application  to  the  demonstration  of  con¬ 
cepts  in  the  previous  sections  and  analysis  of  errors  in  this 
section.  The  study  of  the  theory  of  Poisson  processes  has  also 
provided  the  concepts  necessary  for  digital  simulation  of  the 
LV  signals  and  their  detection  by  photon  counting  systems. 

The  concepts  and  a  FORTRAN  computer  program  which  we  have 
developed  for  this  simulation  are  provided  in  Appendix  B.  The 
simulation  provides  a  method  of  investigating  such  nonlinear 
effects  as  processor  dead  time  and  counter  saturation  which  we 
have  not  yet  been  able  to  do  with  analysis. 

Bias  Errors 

The  sources  of  bias  errors  which  have  been  studied  in 
the  measurement  of  mean  flow  are  as  follows:  fringe  number, 
turbulence  intensity,  time  smear,  and  dead  time. 

Fringe  number  and  turbulence  intensity.-  When  the  number 
of  fringes  N  in  the  probe  volume  is  small,  subtraction1"  of 
R^(x  -  3Tm/4)  -  R^(t  "  Tm/4)  is  negative  instead  of 

zero.  The  error  becomes  significant  only  for  small  values  of 
Nf,  The  signal  frequency  estimate  using  the  zero  criterion 
will  be  too  small.  We  have  computed  the  error  in  percent  in 
the  following  manner. 

Let  R^(3qAx)  -  R^(qAx)  A  ^(qAi).  For  large  number  of 

♦ 

This  error  pertains  to  the  dual  correlate  and  subtract 
technique. 
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T 

fringes  (N**»)  then  iKqAx)  ■  0  when  At  *  ^  where  Tm  is  the 
mean  signal  period.  When  N<*>,  then 

Rx(3Tffl/4)  -  Rx(Tm/4)  -  Ae  -  <J»(Ax  -  Tm/4>  <  0  (61) 

The  correct  value  of  qAx  for  i|>(qAx)  ■  0  is  qAx  ■  Tm/4  +  W r7r"7T) 

^  v  m' 

and  the  fractional  error  e  in  accepting  iKqAx)  -  0  as  indica¬ 
tion  that  qAx  ■  Tm/4  is  approximately 


e 


4MTm/4) 

V'<V4> 


(82) 


We  have  evaluated  e  for  various  levels  of  p  *  rms  turbulence/ 
mean  velocity  and  for  N.  -  number  of  static  optical  fringes  in 

2  i 

1/e  signal  width. 

The  result  of  the  parametric  computations  are  presented 
in  two  forms  in  Figures  6  and  7.  From  the  results  we  conclude 
that  the  bias  error  is  not  very  sensitive  to  turbulence 
intensity  at  the  levels  shown.  A  first  order  correction  inde¬ 
pendent  of  turbulence  is  thus  possible.  A  second  order  correc¬ 
tion  is  possible  with  only  very  imprecise  estimates  of  turbu¬ 
lence  intensity. 

Time  smear  error  and  correction  techniques.-  The  photons 
represented  by  the  count  n^  at  time  kAx  were  in  fact  smeared 
in  occurrence  time  over  the  range  [(k-l/2)Ax,  (k+l/2)Ax], 

The  effect  of  this  is  usually  neglected  by  assuming  Ax  is  small 
compared  with  any  significant  variations  in  the  classical  sig¬ 
nal.  This  is  never  true  in  the  most  difficult  experimental 
cases  where  the  electronic  speed  limitations  force  Ax  to  be  an 
appreciable  function  of  the  signal  period  Tm>  The  result  pro¬ 
vided  in  equation  (24)  for  photon  correlation  included  the  time 
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smear  effect  in  the  triangular  function  A  which  convolves  with 
the  autocorrelation  of  the  classical  signal.  After  equation 
(26)  the  triangular  function  was  treated  as  a  delta  function, 
i.e.,  in  the  limit  of  small  At, 

rf  i(t/iT)  +  j(t)  css) 

and  equation  (27)  results  from  the  convolution. 

Now,  if  we  do  not  make  the  limit  assumption, *  we  note 
that  the  effects  of  the  triangle  function  are  easily  displayed 
in  the  frequency  domain:  From  Bracewell  [24]  we  have  the  rela¬ 
tionship  , 


00 

/  ^  A(T/AT)e’J2lTf  «  sine2  (fAt)  (64) 

00 

where 


sine  (fAT) 


sin  irf At 
irf  At 


(05) 


This  Fourier  transform  relationship  is  illustrated  in  Figure  8. 
The  convolution  theorem  of  Fourier  transform  theory  assures  us 
that  in  the  frequency  domain,  the  effect  of  the  convolution  in 
equation  (24)  is  a  product.  In  other  words,  the  frequency 
spectrum  associated  with  the  signal  correlation  by  Fourier 
transform  is  attenuated  by  a  low-pass  filter  whose  form  is 
sine  (fAT).  This  function  is  plotted  in  Figure  8c.  As  the 
figure  shows,  there  is  little  attenuation  when  ffflAT  -  AT/Tm  ■  1/16, 
i.e.,  when  the  clock  period  is  1/16  of  the  signal  period.  For 

^  . '  .  -  -  -  - .  “jr . " 

The  effect  of  time  discretization  is  a  sine  (f )  type  of 
low-pass  filter.  This  is  similar  to  a  result  in  our  previous 
work  of  determing  turbulence  power  spectra  from  randomly  tried 
samples  [23] . 
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very  low  turbulence  and  a  large  number  of  fringes  Nf,  the 

signal  spectrum  would  be  small  In  width  and  At  could  be  made 

as  large  as  Tffl/2.*  However,  for  values  of  At  >  Tm/8  we  see 

that  a  correction  of  the  velocity  probability  density  function 

could  be  useful  for  all  but  very  small  turbulence  levels. 

Such  a  correction  would  be  effected  by  multiplying  the  trans- 

2 

form  of  the  correlation  estimate  by  1/sinc  (ATf)  prior  to 
final  data  interpretation. 

Dead  time  effects.-  No  physically  realizable  photo¬ 
detector  and  electronic  counter  combination  can  be  constructed 
without  some  dead  time;  i.e.,  a  period  of  time  following  a 
threshold  crossing  by  the  analog  photodetector  voltage  wave¬ 
form  during  which  no  additional  crossing  events  will  be 
counted.  This  dead  time  is  typically  10  nsec  for  commercial 
photon  discriminator  circuits  at  the  time  of  this  writing. 

There  is  no  fundamental  reason  why  this  cannot  be  reduced  to 
less  than  Snsec  with  the  fastest  photomultiplier  tubes  and 
counting  circuits  now  available.  We  will  also  distinguish  two 
other  types  of  counting  dead  time.  The  first  of  these  is 
"pulse  pile  up"  in  which  the  photodetector  analog  waveform 
remains  above  the  threshold  level  due  to  there  being  more  than 
one  photoelectron  event  within  the  pulse  response  time  of  the 
photo  detector.  There  is  also  a  brief  interval  during  the 
periodic  counting  interval  during  which  the  counter  is  being 
reset  and  is  not  available.  This  is  true  even  when  two  counters 
are  used  alternately,  since  a  small  guard  interval  is  required 
to  keep  both  counters  from  counting  the  same  event.  Practical 
dead  intervals  Tr  may  be  on  the  order  of  2nsec  or  less. 


*This  limit  could  be  actually  exceeded  because  the 
sampling  theorem  roally  refers  to  bandwidth,  not  maximum  fre 
quency;  however,  the  attenuation  would  then  be  severe. 
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The  counter  dead  interval  t  will  be  considered  first 

r 

since  it  is  the  simplest.  To  first  order  it  can  be  neglected, 
but  correction  in  the  formulas  is  simple.  It  is  only  neces¬ 
sary  to  replace  At  by  Ar-tr  in  the  theoretical  formulas  where 
the  At  refers  to  the  counting  width.  The  amount  of  delay 
remains  equal  to  some  integer  multiple  of  At.  As  a  simple 
example  the  mean  value  of  nk  becomes  <X(t )>(Ax-Tr)  instead  of 
<X(t )>At . 

The  effects  of  pulse  pile  up  are  least  known  for  the 
photon-limited  signal  cases.  For  higher  photon  rates  the  PMT 
analog  waveform  will  generally  remain  above  the  photon  counting 
threshold  and  so  produce  no  counting  effects.  The  pulse  pile 
up  effects  can  be  studied  with  the  present  signal  simulation 
program  (see  Appendix  B)  when  a  more  complete  processor  simu¬ 
lation  subroutine  is  completed  in  the  future,  At  present  the 
photon  processor  simulation  is  idealized  so  that  discretized 
photoelectron  event  times  are  used  directly  without  synthesis 
of  the  PMT  anode  waveform  and  a  threshold  crossing  circuit 
model . 

The  effects  of  discriminator  dead-time  effects  have  been 
presented  analytically  by  Jakeman*  for  the  case  of  Gaussian 
optical  electric  field  statistics.  We  have  not  yet  extended 
the  theory  to  the  single-particle  LV  signal  situation  (non- 
Gaussian  field  statistics)  and  can  offer  no  improvement  analyti¬ 
cally  here  over  Jakeman's  results.  However,  we  dispute  Jakeman's 
conclusion* that  in  practice  dead  time  errors  are  not  a  serious 
limitation  of  photon  correlation.  The  seriousness  of  the 
effect  is  demonstrated  below  using  the  computer  simulation  pro¬ 
gram  provided  in  Appendix  B. 

We  have  simulated  a  wind  tunnel  instrumentation 
example  with  the  constants  given  in  Table  1.  The  results  are 

*Page  118,  [12] . 
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Theory 

Symbol 

Program 

Name 

Value 

Description 

T 

TTOT 

10” 3sec 

Total  time  limit 

MNTOT 

3x  103 

Photoelectron  limit 

\  Snk 

KMAX 

480 

Number  of  At's  per  burst 

A 

*< 

V 

AFAC 

10^/sec 

Mean  peak  pedastal  rate 

HIGH 

10 

Max  Xj  ■  High*<Xj> 

i 

LOW 

0.1 

Min  Xj  *  Low*<Xj> 

Xj  -  AFAC  (not  random) 

Mean  time  between  bursts 

IR2 

4 

r  i/R 

TB 

5x  10"6 

THEORY 

TRUE 

Logical:  Selects  ideal 
processor 

MAXC 

103 

Maximum  count  of  processor 
(high  to  avoid  limitation 
here) 

A 

a/5 

D 

8  x  10“7sec 

1/e^  half  width  of  bursts 

CON 

3 

D  CON  -  Total  burst  width 

L 

AME 

1.0 

Burst  modulation  index 

FO 

6 . 25  x  108Hz 

Signal  frequency 

At  /At 

ITAU 

4 

Ratio  of  processor  resolu¬ 
tion  to  simulation  resolu¬ 
tion 

■ 

IP 

20 

Maximum  correlation  delay 

WAVE 

TRUE 

Logical:  True-Bursts  present 

xb-° 

CONST 

FALSE 

Background:  True-X^  present 

vAt 

DEAD 

0,1, 2,3,4 

Dead  time  in  At  units 

At 

DT 

5  nsec 

Simulation  resolution 
interval 

At 

DT*ITAU 

20  nsec 

Processor  resolution 
interval 

ONE 

FALSE 

Logical:  Time-Dual  Correlate 

TABLE  1.  PARAMETER  SELECTIONS  FOR  DEAD  TIME  SIMULATION. 
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plotted  in  Figure  0.  For  comparison,  the  figure  also  shows 
the  theoretical  expected  value  of  the  accumulator  sums  as 
predicted  by  the  following  equations.  We  use  the  values  in 
Table  1  with  equations  (40),  (41)  and  (44)  to  obtain  the 
mean  level  <X(t)>  as 

CL 

<X>  ■  2.005  x  10  photoelectrons/sec  (06) 


<Mp>-  80.4  +  283.0  (1  +  \  cos2Trfmr )e  2(B-60x  10  )  (67) 


where 


-9 

t  ■  (pAt)  »  p(20xl0  sec) 

f  -  0.25  x  106  Hz 
m 

From  Figure  0  it  is  dear  that  for  dead  times  which  are  an 
appreciable  fraction  of  the  clock  interval,  there  will  be  dis¬ 
tortion  of  the  correlation  functions.  This  does  not  appear  to 
be  a  problem  for  At  ■  Tm/8  except  at  the  first  delay  valve  and  a 
general  amplitude  reduction  but  further  study  is  needed.  The  dead 
time  effect  is  seen  to  seriously  affect  the  first  delay;  this 
will  seriously  affect  the  dual  correlate  approach  if  At  -  Tm/4. 

Variability  Error 

In  this  section  we  derive  formulas  for  the  fractional  rms 
error  of  a  mean  flow  measurement  with  a  dual  correlate  and  sub¬ 
tract  system  in  terms  of  the  mean  photon  rate  and  the  system 
response  (integration)  time.  We  assume  that  the  accumulator 
sum  is  zero  after  summing  for  total  time  T.  This  implies  that 
the  clock  frequency  is  in  error  by  an  amount  required  to  cancel 
a  random  (finite  time)  variability  error  in  the  accumulator  sum. 
These  compensating  errors  are  assumed  small  so  that  a  perturba¬ 
tion  may  be  used. 


4b 


An  error  in  the  accumulator  sum  is  equivalent  to  an 
error  eT  in  the  estimation  of  one  quarter  of  the  mean  signal 
obtained  by  dividing  by  the  slope  of  the  discriminator  function 
A(t) 


eT  *  VA'<V4> 


(68) 


where  A(t)  is  an  abbreviation  for  <Mpq>  given  earlier  in  equa¬ 
tion  (56)  and  (57).  The  fractional  error  e  in  estimating  the 
period  Tm  is  thus 


e 


eA/A,<y4> 


4ea 


(89) 


This  is  the  same  form  as  equation  (62)  with  the  only  difference 
being  the  type  of  error.  We  proceed  by  obtaining  a  simplified 
expression  for  A(Tm/4)  by  neglecting  the  Gaussian  exponentials 
in  equation  (57),  Evaluation  of  the  derivatives  gives 

A'(Tffl/4)  »  NAt2R<X ?>  a  (70) 

J  m 

The  quantity  is  the  rms  deviation  of  the  accumulator  value 

after  data  collection  time  T.  The  evaluation  of  this 
PQ 

quantity  is  discussed  in  Appendix  C  and  the  result  for  cases 
where  the  steady  background  light  is  larger  than  the  signal  is 
obtained.*  For  this  case  we  have 

eA  -  /"SIT  IpAt  (71) 

When  these  results  are  combined  we  have  the  fractional  error 
in  the  velocity  estimate  given  by 

*This  is  simpler  to  evaluate  than  the  low  background 
case  and  gives  a  bound  on  the  time  required  to  produce  a  given 
error  in  cases  where  the  background  is  less  . 
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e 


2*b 

oi/TAt  R<Xj>tt/u 


(72) 


where 

X^  ■  background  photoelectron  rate  >>  signal  rate 
a  ■  burst  1/e  hall  width 

2 

<X.>  ■  mean  square  peak  pedastal  burst  amplitude 
J  (photoelectron  rate) 

R  »  rate  of  occurrence  of  bursts 

T  »  total  time  ■  NAt 


ADVANCED  PHOTON  PROCESSOR  FOR 
TURBULENCE  TIME  STATISTICS 


Introduction 

In  this  section  we  propose  and  analyse  a  photon -processing 
scheme  for  estimation  of  the  temporal  autocorrelation  of  the 
time  varying  velocity  fluctuations.  Part  of  the  basis  for  this 
is  the  fact  that  the  frequency  discriminator  characteristic  of 
Figure  5  applies  not  only  to  the  long  time  average  of  the  dual 
correlate  sud  subtract  sum  8pq  but  also  Instantaneously  in  a 
conditional  statistical  sense.  We  way  therefore  tune  At  to  the 
value  which  centers  the  long-term  average  at  the  zero  of  the 
frequency  characteristic  and  then  obtain  short  periodically 
occurring  accumulations  of  ft  whose  conditional  expected 
values  follow  the  velocity  deviation.  Even  when  the  photoelec¬ 
tron  rate  is  small,  the  correlation  of  the  velocity  fluctuations 
may  then  be  obtained  by  conventional  digital  correlation  of 
the  sequence  of  short  time  ftp(J’s. 

We  define  here  AT  as  a  period  greater  than  At,  the 

processor  counting  interval,  and  less  than  the  significant 

times  of  change  of  the  turbulent  fluctuations.  The  quantity 

Sl is  a  sum  over  the  interval  [(n-l)AT  <  t  <  nAT]  of  m, 
pqn  k 

defined  in  equation  (55)  as 
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(73) 


®k  "  nk(nk-p  “  nk-q); 

i . e. ,  ft  -  ft _  over  the  nth  AT  interval.  We  will  consider 

pqn  pq 

an  autocorrelation  estimate  of  the  ftpqn  sequence: 

V1 

V14T)  ■  lml  flpqnftpq(nti)  <74> 

We  will  show  that  the  expected  value  of  this  sequence  contains 
the  shape  of  the  turbulence  autocorrelation  function  under 
certain  conditions.  Further  it  appears  that  the  magnitude  of 
the  turbulence  intensity  may  be  obtained  by  a  normalization 
procedure  which  will  be  discussed. 

Spectral  Analysis  of  Randomly  Sampled  Signals 

Continuous  Functions.-  In  our  previous  development  [23] 
we  showed  that  correlations  and  frequency  power  spectra  can 
be  obtained  from  randomly  timed  discrete  samples  of  a  velocity 
component  u(t)  where  the  sampling  function  z(t)  is  a  uniform 
Poisson  impulse  process.  The  same  principles  may  be  general- 
ized  to  the  present  more  complicated  data  processing  problem. 
First,  let  us  assume  that  a  continuous  sigral  s(t)  is  available 
such  that 


s(t)  -  x(t)p(t) 

(75) 

p(t)  -  u(t)/U 

(76) 

where  u(t)  is  a  zero-mean  time-varying  velocity  component 
deviation  from  the  mean  component  0  and  the  random  process 
x(t)  is  a  filtered  Poisson  shot  noise  process  (see  Appendix  A) 
statistically  independent  of  u(t): 


49 


(77) 


X(t)  -  l  b  . f(t  -  T  .  ) 

>  j«co  J  J 

Here  f(t)  is  a  low-pass  impulse  response  function,  the  Tj's 
are  random  occurrence  times  which  obey  a  stationary  Poisson  * 
law,  and  the  bj'B  are  equally-distributed  statistically- 
independent  random  amplitude  variables.  The  functions  are 
illustrated  in  Figure  10.  From  the  properties  of  Poisson 
processes  we  know  that  if  f(t)  is  a  positive  function,  then 
R  (t)  ■  <x(t)x(t + t )>  is  a  positive  function  and,  making  use 
of  the  independence  assumption 

Rp(T)  -  <p(t)p(t  +  t)>  (78) 

-  Rs(t)/Rx(t) 

We  observe  that  the  zero  value  of  Rp(T)  is  the  normalized 
mean-square  turbulence  intensity  <p2>.  Figure  10  illustrates 
the  fact  that  when  the  duration  of  the  pedastal  correlation  is 
small  compared  with  the  duration  of  the  velocity  correlation, 
the  value  <p  >  may  be  obtained  approximately  from  Rp(r)  where  t 
is  small  but  greater  than  the  pedastal  duration.  If  the  turbu¬ 
lence  intensity  is  obtained  in  some  other  manner,  we  may  theo¬ 
retically  obtain  the  shape  of  R  (t)  even  without  obtaining 
R  (x),  except  in  the  vicinity  of  t  ■  0,  since  R  (t)  -  constant 

X  X  . 

elsewhere . 

We  may  extend  the  above  reasoning  to  the  situation  where 
two  sets  of  processes  are  available 

sx(t)  -  x^tJpjCrj  ,t)  (70) 

s2(t)  "  x2^t^p2^2’t^ 

So  long  as  x^  and  Xg  are  independent  of  p^  and  p2,  then 
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Rs12(t )  "  <s1(t)s2(t+T)> 


(80) 


-  <x1(t)x2(t+T)xp1(r1,t)p('r2,t+T)> 

*  Rx12<t>Rp12(t) 

and  the  cross  correlation  may  be  obtained  by  division 

WT>  '  Rsl2<T>/Rxl2'T>  <81> 

The  frequency  power  spectra  are  Fourier  transforms  of  the 
correlation  functions. 

Sampled  Data  Estimates.-  Zn  the  preceding,  we  assumed 
s(t)  and  x‘(t)  were  continuous  and  computed  statistical  corre¬ 
lation  functions.  We  now  suppose  that  real-valued  discrete 

SB  and  Xn  of  .(t)  and  x(t)  are  taken  at  uniform  inter 
vale  AT  which  are  small  compared  with  the  duration  of  the 
random  sampling  pulse  f(t)  and  correlation  time  of  p(t).  We 
may  form  discrete  estimates  of  R_(t)  and  Rv(t)  at  t  ■  iAT  by 
computing  a  finite  time  average,  for  example 

Nt-i 

j,  8nSn+l  '82> 

where  N^AT  is  the  accumulation  time  T  and  the  expected  value 
ie 


<Rs(iAT)>  -  Rs(x  -  iAT)  (83) 

If  such  estimates  are  formed  for  R  and  R  ,  then  we  may  form 

A  o 

an  estimate  for  Rp(x)  by  division: 


Rp( iAT)  £  Ra(iAT)/Rx(iAT) 


(84) 


This  procedure  does  not  insure  that  <Rp(iAT)>  is  an  unbiased 
estimate  but  it  becomes  one  in  the  limit  of  sufficiently  long 
accumulation  time  when  the  numerator  and  denominator  converge 
to  their  respective  expected  values. 

Conditional  Expectation  Again.-  The  problem  is  still 
more  complicated.  Because  we  do  not  have  the  classical  optical 
signal  X(t)  available  for  direct  observation,  we  cannot  form 
x(t)  and  s(t),  or  even  Xn  and  Sn>  We  will  process  the  fast 
phot on- counting  sequence  nk  to  obtain  discrete- valued  sequences 
Mxn  and  Msn  whose  conditional  expected  values  given  X(t)  are 
proportional  to  X„  and  Sn. 

We  thus  have  an  estimate  fig(iAT)  defined  not  by  equation 
(82),  but  by 

N t -i 

CaB,(lAT)  -  nZi  MgnMs(n+1)  (88) 


The  expected  value  of  the  estimate  is 


<UBnMs(n+l)>  ■  «K«jMt>><M,(n+1)|XCt)»  (86) 


i  C%3  3  > 

c  J  J+i 


-  C  R  (iAT) 

s 


where  C  is  the  proportionality  constant.  This  result  depends 
on  the  fact  that  given  X(t)  the  discrete  random  variables 
Msn  and  Ms(n+i)  formed  from  non  overlapping  portions  of  the 
sequence  nfc  are  conditionally  independent.  This  follows  from 
the  fact  that  n^  is  Poisson  counting  process. 


In  equation  (86)  the  notation  <x|X>  denotes "expected 
value  of  x  given  X"  and  we  made  use  of  the  fact  that  the  uncon¬ 
ditional  expectation  may  be  obtained  in  two  steps:  first,  with 
respect  to  x  given  X;  then  with  respect  to  X. 

Application  to  Photon  Processing  of 
Turbulence  Correlation 

It  remains  for  us  to  identify  measurable  quantities  with 
the  properties  attributed  to  Msn  and  Mxn  above  and  to  describe 
a  means  of  implementation. 

Conditional  frequency  discriminator.-  We  now  reconsider 
the  estimator  Mp(l  defined  in  equations  (54)  and  (55).  We 
restrict  ourselves  to  the  signal  model  given  in  equations  (30) 
and  (38)  with*  Xfe  ■  0  and  with  rarely  overlapping  bursts; 
i.e,,  a<<l/R.  This  is  the  "low  density"  shot  noise  case  dis¬ 
cussed  by  Papoulis  tlOl,  page  574.  For  "low  density"  shot 
noise  we  may  use  the  approximation  that  [10]  if 


i < t )  ■  Zh(t  -  tj) 


then 


s2(t)  *  Zh2(t  -  Tj) 
with  these  constraints  we  obtain 

<hknk.pU<t)>  -  Af2X(t )X(t  -  p Ax) 

*  At2  l  Xj2f2(t  -  t j )  <1  +  cos  Wj(t-Tj) 

J“— 06 

+  cos  Wj(t  -  Tj  -  pAt)  +  cos(2ojjt  -  2u)jT j  -  pwjAr) 
+  cos  pwjAt} 


(87) 

(88) 


(80) 


the 


The  Xw 
subtract!. 


■  0  assumption  does  not  affect  the  result  due  to 
on . 
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In  equation  (80)  the  product  of  cosine  terms  has  been  expanded 

as  sum  and  difference  frequency  terms,  and  the  delay  pAt  has 

been  neglected  in  the  envelope  function  f(t).  With  mk  defined 

by  equation  (55)  we  may  determine  <mk|X(t)>  from  two  sets  of 

terms  like  that  in  (80).  We  then  determine  <KL  >  as 

pqn 


a  nL 

<MDanU(t)>  "  I  <mJX(t)> 
pqn  k-(n-l)L  K 

nAT 


(90) 


nnT  t  » 

T  J(n-l)iT  ( “  TJ > tcOB  P“jAT 
-  cos  quij At]  +  [Band  pass  termslj  dt 


where 


L  -  AT/At 


(91) 


In  the  right  hand  side  of  this  equation  the  sum  is  replaced  by 
a  time  Integral  which  it  approximates.  Since  AT  is  long  com¬ 
pared  with  the  Doppler  period,  the  bandpass  terms  average  to  a 
small  value;  however,  since  AT  was  assumed  short  compared  with 
the  burst  envelope  duration,  the  integral  does  not  smooth  the 
first  expression  in  the  right  hand  side  of  (90).  The  result 
is  therefore 


2  2 

Xj  f  (t  -  Tj ) [cos  pATWj  -  cos 


qAiw j]  (02) 


m  .4^4?  X2Jlp(t)[cos  p A t u ( t )  -  cos  qATw(t)] 

By  proper  selection  of  a  high  speed  digital  clock  period, 
At,  the  sequence  approximates  the  sequence  &sn  described 
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above  for  certain  values  of  the  integers  p  and  q.  For  other 
choices  of  p  and  q  the  sequence  $  approximates  the  sequence 
llxn<  Thus,  velocity  correlations  may  be  obtained  by  simul¬ 
taneously  duplicating  the  computations  of  (73)  with  two  differ¬ 
ent  sets  of  integers  p  and  q,  and  then  proceeding  with  soft¬ 
ware  processing  of  the  two  sequences  MBn  and  MXR.  To  show 
this  last  link  in  the  procedure  we  examine  the  cosine  differ¬ 
ence  term  in  (92).  For  the  MflQsequence  we  select  p  and  q  and 
adjust  the  variable  At  to  the  conditions  giver,  previously  for 
the  mean  frequency  discriminator  in  equations  (58)  and  (59). 

(P  -  <l)wmAt  -  ft  (93) 

« 

pu)mAx  ■  3ft/2 
qu>mAT  -  ft/ 2 

Example  selections  for  p,  q  and  At  are  3,l,Tm/4;  0,2,Ta/8; 
12,4,Tm/16;  where  Tm  ■  2tr/wm  is  inverse  of  the  mean  signal 
frequency.  Expansion  of  the  bracketed  cosine  term  in  (92) 

gives 

HftAojj  3ftAu)j  *| 

2^-J-)  +  sin(^j — “-)J  (94) 

where  we  have  used  the  substitutions 


Equation  (94)  was  shown  plotted  in  Figure  5  since  it  is  also 
the  approximate  long-time  average  frequency  transfer  charac¬ 
teristic  for  the  mean  flow  measurement.  The  trigonometric 
substitutions  along  with  the  small  angle  approximation  sin  0-0 
show  that  for  small  turbulence  levels,  equation  (92)  becomes 

k 

I  X(t  )>  “  ATAtX2  (t)irp(t)  (96) 

Where  p(t)  *  durini  non-zero  portions  of  XAp(t). 

Normclizing  sequence.-  Under  the  same  selection  of  At  as 
in  (93)  but  with 

pu>mAT  *  2tt  (97) 


qwmAt  ■  tt 


a  different  effect  is  obtained.  For  example,  if  Man  is  obtained 
with  p  «■  3,  q*l  and  we  let  p*4,  q  ■  2  in  computing  Mxn,  the 
result  corresponding  to  (94)  is 


[cos  pATWj  -  cos  qAtWj] 


(98) 


The  small-angle  approximation  for  the  cosine  function  is  unity, 
so  the  result  corresponding  to  (96)  is 

<MxnU(t)>  “  ATAT\2Jlp(t)  (99) 

When  we  note  that  Auj/wm  ■  u j / 1?  and  that  the  ATAt  term  cancels 
by  division,  we  see  from  equations  (96),  (99)  and  (86)  that 
the  autocorrelation  of  M  divided  by  the  autocorrelation  of 
produces  approximately  it  <u(t  )u(t  +  t)>/ TJ2 .  Except  for  a  fac¬ 
tor  of  C2  *  tt2,  the  result  is  normalized  in  such  a  manner  that 
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fractional  mean-square  turbulence  intensity  is  directly  obtain¬ 
able  from  the  first  lag  value  of  the  final  autocorrelation 
estimate  or  a  projection  back  to  the  zero  delay  value.  (The 
analysis  we  have  presented  is  not  valid  for  the  zero  delay 
value  without  separate  consideration.  It  is  possible  that  cer¬ 
tain  effeots  will  canoel  and  cause  this  to  be  a  valid  point 
also.) 

Practical  considerations.-  The  approach  outlined  above 
for  normalization  would  require  as  much  hardware  to  measure  the 
Mxn  **4uence  as  the  Mfln  sequence.  In  addition,  the  small  angle 
approximation  whioh  was  applied  to  (98)  is  not  valid  over  a 
very  large  range  of  velocity  deviation.  Ve  have  conceived  two 
other  less  expensive  approaches.  The  first  of  these  consists 
of  eliminating  the  second  Mxn  channel  completely  and  relying  on 
the  fact  that  the  shape  of  the  correlation  Rp(T)  is  approximated 
by  correlating  M§n  except  in  the  vicinity  of  a  burst  duration 
from  the  delay  origin.  If  it  were  desirable  to  normalize  the 
function  at  the  rest  of  the  delay  locations,  we  could  evaluate 
the  required  division  oonstant  as 

ATaAtair2  <X2Ap(t)>2  -  <Rg(nAT)>,  nAT  »  a  (100) 

o 

where  we  have  already  discussed  a  possible  estimate  for  <a  (t )> 
in  equation  (46).  Either  this  estimator  or  the  M  approach 

AM 

described  above  could  be  computed  with  a  single  channel  elec¬ 
tronic  system  by  electronically  switching  the  delay  variables 
p  and  q  and  storing  the  {ftpqn)  sequences  in  different  portions 
of  a  computer  memory  before  computing  the  second  level  corre¬ 
lations.  Time  sharing  like  this  would  not  be  as  efficient  in 
possible  cancellation  of  some  of  the  statistical  variability 
error,  but  significant  cost  savings  would  result. 
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DISCUSSION 


System  Design 

We  have  provided  a  specification  for  an  advanced  computer 
controlled  photon  counting  processor  in  Appendix  D.  The  design 
allows  the  system  to  be  operated  as  a  dual-correlate  and  sub¬ 
tract  processor  or  sequentially  as  a  correlator  by  time  sharing 
the  multiplier.  The  system  provides  synchronous  3-bit  X  3-bit 
operations  at  up  to  100  MHz.  The  design  uses  "slow"  emitter 
coupled  logic  (ECL)  which  is  optimum  for  the  speed  range  speci¬ 
fied.  The  system  may  be  operated  as  an  advanced  processor  by 
dumping  the  accumulator  values  to  computer  memory  at  rates 
limited  only  by  the  computer.  The  design  allows  two  identical 
units  to  be  used  together  for  either  velocity  cross-correlation 
measurements  or  for  simultaneous  second-channel  normalization. 

An  analog  feedback  loop  is  provided  for  automatic  mean  velocity 
acquisition.  In  addition,  the  system  can  bo  used  to  measure  any 
type  of  multiple-interval  photon  statistics  by  selecting  the 
correlator  clock  period  (continuously  variable  both  manually  and 
electronically)  and  using  the  computer  to  sample  the  values  of 
nk  stored  in  the  delay  line. 

Sensitivity  Comparison 

There  are  four  primary  sources  of  variability  error. 

These  are  the  random  turbulent  flow  itself,  the  random  occur¬ 
rence  times  of  the  scattering  particles,  the  random  scattering 
cross  sections,  and  the  random  time  of  photon  events.  Photon 
correlation  methods  are  linear  in  the  sense  that  the  effects  of 
two  simultaneously  occurring  scatterers  are  added.  This  is 
beneficial  in  that  it  avoids  the  non-linear  zero-crossing  cap¬ 
ture  effect  of  classical  FM  systems  and  thus  the  error  problems 
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of  multiple  scatters;  it  is  not  beneficial  in  the  sense  that 
random  amplitude  effects  will  contribute  to  variability  error, 
but  these  distinctions  have  not  been  carefully  analyzed.  The 
error  versus  data  collection  time  trade-offs  which  result  due 
to  the  random  occurrence  times  of  Ecatterers  and  due  to  mini¬ 
mum  limits  due  to  the  turbulence  itself  have  been  previously 
analyzed  for  burst  counters.  That  result’1'  may  be  expressed  as 

T "  ^ +  <101) 

■  i 

where  B  is  the  equivalent  power  bandwidth  of  the  turbulence, 

2 

R'  is  the  mean  rate  of  accepted  signal  bursts,  and  <p  >  and 
<e  >  are  the  normalized  mean-square  turbulence  intensity  and 
mean-square  estimate  error,  respectively: 

<p2>  .  (102) 

Tr 

<e2>  . 

this  result  indicates,  for  example,  that  for  a  turbulence  equiva¬ 
lent  power  bandwidth  of  2KHz,  10%  turbulence  intensity,  and  0.3% 
desired  rms  error,  the  data  collection  time  would  be  at  least 
0.28  seconds  even  if  the  continuous  signal  U(t)  were  available 
for  processing,  and  would  be  considerably  longer  if  R'  were 
less  than  4000/sec. 

In  the  following  we  examine  the  implications  of  equation 
(72)  for  transonic  wind  tunnel  measurements.  In  order  to  include 
the  mean  value  of  the  signal  bursts  along  with  the  assumed 

j|l 

page  10,  equation  (2.14),  reference  123]. 
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background  light  in  equation  (72)  we  use  equation  (40)  for 
<X(t)>  and  (72)  becomes 


e 


« 


»<vx.> 

avTZ?  R <X j >  7T /F 


where  the  mean  rate  due  to  signal  bursts  is 


(103) 


Xs  "  R<Xj>  ^  0  (104) 

which  is  less  than  the  mean  peak  rate  <Xj>  when  R  is  less  than 

the  inverse  of  the  effective  burst  width,  /tr  ot,  Now  by  defin- 

2 

ing  the  variance  cr^  of  the  pedastal  peak  rates  Xj  as 

2  2  2 

V  "  <XJ>  -  v 

we  may  rearrange  (103)  as 


(105) 


Although  the  parametric  behavior  of  equation  (105)  is 

intuitively  acceptable  in  other  ways,  the  presence  of  the  term 
2  2 

[1  +  a .  / <Xj >  ]  in  the  denominator  seems  a  little  strange. 
Increased  variability  of  the  scattering  cross  sections  of  the 
particles  would  not  intuitively  decrease  the  mean-flow  varia¬ 
bility  error.  It  is  possible  that  this  indicated  behavior  is 
a  consequence  of  ignoring  the  variability  of  the  classical  power 
X(t)  in  equation  (71)  and  has  no  physical  reality.  However, 
there  does  not  appear  to  be  anything  wrong  with  the  derivation 
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for  the  limiting  case  of  high  background  light  fo  the  strange 

result  may  be  true.  In  any  event,  a  conservative  bound  is 

obtained  by  removing  the  bracketed  expression  in  the  denomina- 

2 

tor  to  obtain  an  expression  not  requiring  knowledge  of  : 


E 


s 


it<Xj>  /TAt 


(100) 


If  we  evaluate  equation  (106)  with  the  following  assumed  typical 
conditions  for  a  transonic  wind  tunnel  measurement  we  obtain 
a  required  data  collection  time  of  0.5  seconds  for  a  1%  rms 
error : 


Xb  -  107  (107) 

Xe  -  106 
<Xj>  -  107 
At  -  10*'8 

In  this  example,  the  mean  signal  photoelectron  rate  is  ten 
times  loss  than  the  mean  background  photo-electron  rate  and 
is  equal  to  the  average  peak  envelope  rate.  With  this  much 
background  light,  the  assumption  of  constant  X(t)  should  be 
valid  with  respect  to  the  photon-fluctuation  induced  varia- 

7 

bility.  The  selection  of  mean  peak  rate  at  10  means  that 

Q 

occurrences  of  photoelectron  count  rates  greater  than  10  /sec 
(the  limit  of  current  hardware  state  of  the  art)  will  be  rare 
and  the  effects  of  nonlinearity  negligible).  The  selected 
ratio  of  <^j>/^s  “  10  implies  that  the  measurement  volume  is 
only  assumed  to  contain  a  scatterer  1/10  of  the  time  on  the 
average. 
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Even  though  dramatic  improvements  over  the  previous  exam¬ 
ple  situation  may  result  from  reduced  background  light,  the 
result  still  indicates  that  practical  measurements  may  be 
obtained  with  a  peak  signal  photoelectron  rate  of  107  sec.  In 
order  to  compare  this  with  the  performance  of  a  burst-counter 
system,  we  must  assume  values  for  p,  q  and  fm,  the  Doppler  fre¬ 
quency.  With  p  *  3,  and  q  ■  1  the  Doppler  frequency  fm  is 
1/4At  ■  25  MHz.  This  would  result  from  If  ■  304 . 8  m/sec  with 
an  optical  sensitivity  of  82  kHz/m/sec.  This  peak  photoelec¬ 
tron  rate  assumed  is  thus  0.4  photoelectrons/cycle  in  the  pres¬ 
ence  of  0.4  background  photoelectrons/cycle.  For  comparison  we 
note  that  standard  optical  noise  formulas,  given  X(t)  » 

Xj(l  +  cos  u)mt)  +  X^,  would  result  in  peak  SNR  of 


X 

SNR  -  \h)  <l08> 

In  our  example,  B  -  25  MHz  and  Xb  -  107;  if  we  choose  Xj  -  109, 
a  factor  of  100  greater  than  in  our  example,  then  the  SNR  is 
10  at  the  peak  of  the  signal  burst  and  1.35  at  the  1/e  signal 
envelope  points.  Since  this  example  represents  marginal  or 
inadequate  SNR  for  burst  counter  operation  we  deduce  that,  even 
with  100  times  more  scattered  power,  only  the  larger-than-aver- 
age  scatterers  would  contribute. 

Under  conditions  of  less  background  light,  the  burst- 
counter  analysis  would  not  be  improved;  however,  the  photon 
counting  system  results  are  expected  to  improve  considerably. 
Thus  we  conclude  that  mean-flow  measurements  with  from  100-1000 
times  less  optical  power  are  feasible  with  the  photon  counting 
system. 
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Future  Work 


During  this  contract  effort  we  have  developed  simulation 
techniques  which  can  be  extended  to  be  applicable  to  all  types 
of  LV  signal  processors  for  any  level  signals.  The  program  in 
Appendix  B  only  represents  the  beginning  of  how  these  techniques 
can  be  exploited,  and  we  did  not  have  time  under  the  present 
contract  to  use  that  program  except  for  check-out  waveform  simu¬ 
lations  and  the  dead-time  example.  Similarly,  the  higher  order 
moment  equations  developed  in  Appendix  A  have  not  been  yet  used 
to  extend  the  variability  error  analysis  to  include  the  low 
background  case. 

A  system  such  as  that  specified  in  Appendix  D  should  be 
constructed  and  tested.  The  results  in  this  report  indicate 
that  such  a  system  would  a  low  LV  measurements  to  be  performed 
which  are  not  now  feasible.  In  addition,  the  system  would  be 
a  valuable  research  tool  for  many  other  fields  of  research 
which  require  high-speed  digital  correlation  or  measurement  of 
photon  Interval  statistics. 
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CONCLUSIONS 


The  most  significant  new  results  which  this  report  provides 
are  itemized  below: 


e  A  general  non-st at ionary  Poisson  process  model  for 
dual  scatter  laser  velocimeter  (LV)  signals  and  noise 
valid  from  high  level  signals  down  through  photon 
resolved  signals. 

e  Computer  simulation  algorithms  valid  over  the  entire 
range  of  signal  levels,  which  may  be  used  to  evaluate 
any  new  type  of  LV  signal  processor. 

e  A  description  and  statistical  analysis  of  both  con¬ 
ventional  photon  correlation  and  Dual  Correlate  and 
Subtract  frequency  discriminator  technique  for  mean, 
turbulence  intensity  and  turbulence  correlation 
estimations  from  photon  resolved  signals. 

e  A  system  design  for  an  advanced  photon-counting 
processor  which  implements  both  conventional  photon 
correlation  (sequentially)  and  the  Dual  Correlate 
concepts  with  time  resolution  to  10  nsec. 
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APPENDIX  A 


FILTERED  INHOMOGENEOUS  POISSON  PROCESSES 

Is  this  appendix  we  have  used  the  material  from  Snyder's 
book  [9]  to  derive  the  higher  order  moments  of  a  filtered 
inhomogeneous  Poisson  process  up  through  order  4.  We  begin 
with  introductory  material  similar  to  that  given  by  Papoulls 
[10] . 


Inhomogeneous  Impulse  Processes 

The  input  to  a  random  linear  system  is  an  inhomogeneous 
Poisson  impulse  process  as(t)  given  by 


a(t)  -  l  6(t  -  xt) 


<A1) 


—  00 


where  is  the  set  of  random  occurrence  times,  X ( t)  is  the 
instantaneous  statistical  mean  value  of  z(t),  (and  alBO  the 
inean  rate  of  occurrence  of  the  t^'s),  and  6(t)  is  the  dlrac 
delta  function.  The  random  variables  t1  are  independent  of 
eaoh  other  statistically  and  obey  the  inhomogeneous  counting 
law,  l.e.,  the  probability  of  n»k  occurrences  in  the  interval 
(t1,t2)  is 


P(n(t1#ta)  ■  k) 


(A2) 


where 


U  ""  /  X(t)dt 
t. 


(A3) 


The  quantity  y  is  also  the  mean  and  variance  of  the  random  var 
lable  n(t1,t2). 
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The  Response  of  a  Random  Linear  System-Campbell' s  Theorem 

The  output  of  the  random  linear  system  s(t)  is  the  super¬ 
position  of  the  response  h(t  -  t^,  Y^)  to  each  input  impulse: 

00 

«(t)  -  I  h(t  -  t4,  Yt)  <A4) 

00 

where  {Y^}  is  a  set  of  identically  distributed  independent 
vector  random  variables.  The  random  variable  ?  affects  the 
shape  and  amplitude  of  the  response  function  h(t,7).  In  the 
case  of  the  PMT  signal  it  may  take  the  form  of  a  single  scalar 
amplitude  variable.  In  the  case  of  the  classical  optical  sig¬ 
nal  from  turbulent  flow  both  a  random  amplitude  parameter  and 
one  or  more  random  shape  parameters  due  to  velocity  magnitude, 
direction,  and  probe  volume  translational  entrance  location  may 
be  required.  The  theory  should  be  applicable  so  long  as  the 
set  of  multidimensional  random  variables  Y,^  is  independent  of 
the  set  of  occurrence  times  {t^}.  The  generalized  Campbell's 
theorem  results  for  the  Instantaneous  statistical  mean,  variance, 
and  auto-covariance  of  s(t)  are  given  below,  they  apply  regard¬ 
less  of  whether  individual  pulses  are  resolved  or  not. 

00 

<B(t)>  ■  /  X(t )<h(t  -  t ,Y)>dx  (A5) 

—  00 

oo 

o2(t )  ■  <s2(t )>  -  <s(t )>2  ■  /  X(t ) <h2(t  -  t , Y)>dx  (A6) 

00 

cov[s(t1)s(t2)]  *  <s(t1)s(t2)>  —  <s(t1)xs(t2)>  (A7 ) 

00 

•  /  X(T)<h(t1  -  T,Y)h(t2  -  T,Y)>dT 

00 

where  <>  denotes  expectation  with  respect  to  Y  inside  the 
integral  signs.  For  a  causal  signal  such  as  that  from  the  PMT, 
where  h(t)  is  zero  for  t  >  0,  then  the  upper  limits  of 
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integration  may  be  replaced  by  t  or  the  minimum  of  tj  and  t2 
in  equation  (A7).  For  the  transient  case  where  the  impulse 
signal  z(t)  is  applied  at  t  ■  0,  the  lower  limits  of  integra¬ 
tion  may  be  replaced  by  0. 

Higher  Order  Moments 

Nummary  of  statistics.-  Given  a  filtered,  inhomogeneous, 
compound  Poisson  process: 

00 

s(t)  -  l  h(t,Ti;Yi)  (A8) 

—  00 

where  x^ane  random  occurrence  times  which  occur  with  intensity 
X(t),  and  where  7^  are  vector  random  variables  which  are  sta¬ 
tistically  Independent  and  identically  distributed,  we  obtain 
the  following  result.  The  cumulants  are: 

Yl(ti)  -  r)x  -  /  X(t > <h<t1 , t ;Y>  >  dt  (A9) 

00 

00 

Y12(tl’t2)  "  y12  “  /  ^(T)<i(t1,T;Y)h(t2,T;7)>dT 

00 

00 

Y123(tl,t2't3)  *  /  X(t) <h(t1,T;Y)h(t2,i ;7)h(t3, t ;7)>dr 

00  » 

00 

Y1234(tl,t2,t3’t4)  "  /  X(T)<h(t1,T|7)h(t2,T;7) 

00 

•  h(tg,T;7)h(t4,T;7)>  dr 

The  formulas  which  relate  the  cumulants  to  the  moments 
are  as  follows: 


1.  <e(t1)>  " 


2. 


<s(t1)s(t2)>  -  n12  +  n1n2 


where  y12 


covariance  ■ 


(A10) 
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3.  <a(t1)a(t2)*(t3)>  -  y123  +  iv23  +  n2y13 

1  +  n3M12  +  nln2n3 

4.  <s(t1)e(t2)s(t3)s(t4)>  -  y1234  +  n1n2n3n4 

+  nlY234  + n2Y134  f  n3Y124  +  n4Y123 
+  M12y34  +  p13w24  +  u14w23 
+  *14*2*3  +  *13*2*4  +  *12*3*4 

...  ,V  +  w24*i*a +’ *23ni*4  +  *34*1*2 

>’  1  I  ; 

The  derivation  of  the  above  formulae  follows. 

Derivation.-  The  derivation  of  the  preceding  formulas  is 
straightforward  but  tedious  if  we  are  given  the  joint  charac¬ 
teristic  function 

4uj;i 

+<u)  *  <e^  >  (All) 

where 


(1)  -  (U^,  U>2  t  w3>  w4) 


and 


a  -  ls(t1),s(t2)ls(tg),s(t4)] 


From  Papoulis  we  know  that 


(A12) 


(A13) 
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( A14 ) 


<s(t1)s(t2s(t3)s(t4)> 


'l  d4h(u) 


( j )  du>1 3^2  3Wg  da>4 


a»"0 


and  similarly  for  other  momenta;  i.e.,  we  may  obtain  the  moments 
by  determining  the  appropriate  partial  derivatives  of  the  joint 
characteristic  function  at  u>  «  0  which  is  in  turn  a  vector 
valued  Fourier  transform  of  the  joint  probability  function  for 
the  random  variables  sCt^),  s(t2),  etc. 

The  theory  of  filtered  Poisson  processes  provides  us  with 
the  joint  characteristic  function  of  the  second  kind,  \ p,  where 

¥(w)  -  t,n  <j>  (w)  (A18) 

or 

,  .  ,  :  .  i  i  * 

*<w)  -  (Alfl) 

We  may  therefore  either  write  out  $(w)  and  evaluate  the 
derivatives  directly  or  usfe  the  above  equations  to  first 
express  the  partial  derivatives  of  <J>  (to)  in  terms  of  the  par¬ 
tial  derivatives  of  \Ha>) .  We  have  taken  the  latter  approach 
using  the  chain  rule.  As  an  example 

■  *'  (luJ,Suj'  +  ISj1  <Al7) 

The  rest  of  the  derivatives  are  omitted  here  since  they  get 
progressively  more  lengthy. 

We  now  need  only  to  evaluate  the  products  of  partial 
derivatives  We  have  obtained  at  w  “  0.  From  the  material  given 
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by  Snyder 


♦ 


we  have 


♦<«) 


/  X(x)<eJe(w)  -  1>  dx 


(A18) 


where 

0(oJ)  -  w1h(t1,x;Y)  +  u)2h(t2,x;?)  +  . . .+  wnh<tn,x;Y)  (A19) 
From  the  form  of  the  expression  for  iK<*>)  it  is  clear  that 

?<«)|  *  /  X(x)<e*^  -  1>  dx  *  0  (A20) 

I—  .  -eo 

W"0 

therefore 


4><w)|  -  1 

3-o 


(A21) 


It  is  not  difficult  to  see  that  the  form  of  the  partial 
derivatives  is 


n.i./r: 


3 

tio 


L  ■  Jn  /  *<t)  <h(t±,x;Y)h(t2,x;Y).  .e 

1  2  ...  1  o)-0  -* 


J*<»> 


>dt 


u>»0 

(A22) 


i.e.,  the  partlals  of  i|)(oi)  at  id  ■  0  are  equal  to  the  same  order 
cumulant  except  for  the  jn  factor.  With  this  formula  we  can 
now  go  directly  from  the  expression  of  the  partials  of  <t>  in 
terms  of  the  partials  of  i|>  to  the  desired  higher  order  moments. 
The  result  is  that  given  in  the  summary.  If  moments  higher 

*  '  "  '  . "  ~~  ' 

In  Snyder's  book,  Random  Point  Processes,  equation  (4.15) 
the  lower  limit  of  integration  is  to  correspond  to  the  begin¬ 
ning  of  the  process.  The  upper  limit  is  the  minimum  of  the 
times  t^.  This  assumes  that  h(t1,x;T)  is  a  causal  function. 
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than  the  4th  are  needed  we  would  simply  apply  the  chain  rule 
to  determine  the  higher  derivative  of  <j>  in  terms  of  those  of 

Moments  of  a  Gaussian  Random  Process."  We  have  not  yet 
expanded  all  of  the  moments  of  a  non-zero  mean  Gaussian 
process  for  comparison,  We  observe,  however,  that  the  factori¬ 
zation  property  of  zero-mean  Gaussian  processes  does  not  hold 
for  filtered  Poisson  process.  Otherwise  there  are  great  simi¬ 
larities  except  for  the  cumulants  as 

Gaussian  Poisson  (A23) 

<xx>  ■  <x1>  ■ 

<xlx2>  "  “12  +  nln2  <xlV  "  ^12  +  "ln2 

for  (ni  -  0) 


<x1x2x3>  -  0 

<xlX2x3x4>  -  “12M34  +  ^13^24 
+  ^14^23 


<X1X2X3>  ‘  ^123 
<*1*2x3*4>  -  W12V34  +  V13Va4 
+  U14lJ23  +  T1234 


U-Ol 

4  4 

<x  >  ■  3o 


+ 


0a2n2  + 


<x 


4>  *  3a4  +  6n2a2  +  n4 


+  4n  /  Ui3(t)dt 
+  J\ h4(t)dt 
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APPENDIX  B 


DIGITAL  SIMULATION  OF  LOW  LEVEL  DUAL  SCATTER 
LV  SIGNALS  AND  IDEALIZED  PHOTON  PROCESSORS 

This  appendix  provides  theory  and  an  example  FORTRAN 
program  for  digital  LV  signal  simulation.  Background  theory 
and  several  more  complicated  Poisson  Impulse  simulation  algo¬ 
rithms  are  discussed  by  R*  H.  Forrester,  Jr.  In  a  Masters  Thesis 
[26]  performed  under  Don  Snyder.  The  simple  approach  taken 
here  Is  to  discretize  the  possible  occurrence  times  of  photon 
and  classical  burst  signal  events  to  uniformly  spaced  Inter¬ 
vals  of  length  At,  where  At  is  less  than  any  significant  sys¬ 
tem  integration  time.  This  discretization  of  time  Imposes 
itself  upon  all  temporal  system  parameters,  such  as  processor 
olook  interval  At  and  dead  time  td,  which  would  be  continuous 
variables  in  the  real  situation. 

Theory 

Realization  of  Inhomogeneous  Poisson  impulse  processes.-  A 
sample  function  of  an  Inhomogeneous  Poisson  impulse  process  is 
speolfled  by  a  set  of  event  occurrence  times  (t^)  as  described 
in  Appendix  A.  The  key  procedure  required  for  simulation  of 
LV  signals  Is  therefore  the  generation  of  a  set  (t1>  given  a 
specified  rate  intensity  function  X(t).  This  may  be  done  by 
first  generating  a  homogeneous  (stationary)  process  with  unit 
intensity  X(t)  ■  1  and  then  mapping  the  realization  times 
through  solution  of  an  integral  equation.  Figure  B1  illus¬ 
trates  the  required  mapping  in  a  manner  which  helps  provide 
intuitive  grasp  for  what  follows. 

Let  us  define  the  set  of  interarrival  intervals  {w^} 
between  event  occurrence  times  {t^}  as 
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unit  intensity 


w. 


-  T 


i-1 


(Bl) 


For  &  homogeneous  Poisson  impulse  process  with  constant  rate 
intensity  X(t)-X,  it  is  necessary  and  sufficient  that  {w^ 
be  a  set  of  statistically  independent,  identically  distributed 
exponential  random  variables  with  common  probability  density 
Pw(w),  given  by , 

pw(w)  -  Xe“Xw  (B2 ) 

We  may  therefore  generate  the  set  '  {t^}  of  occurrence  times  for 
a  realization  of  a  constant  intensity  process  as 


t 


i 


Ti-1  + 


w. 


(B3 ) 


4* 

Commonly  available  subroutines  generate  statistically 
independent  realizations  x^  of  a  random  variable  x  uniformly  1 
distributed  on  the  interval  (0,1).  The  exponential  density 
function  is  monotonic  with  an  inverse  function  which  is  commonly 
available  (namely  the  natural  logarithm  function).  It  1b 
therefore  straightforward  to  determine  a  transformation  which 
maps  the  realizations  of  the  unit  uniformly  distributed  varia¬ 
ble  to  the  desired  exponentially  distributed  variable: 


*l  "  T  ln<3T[> 


(B4 ) 


In  summary,  generation  of  unit  uniform  random  variables  x; 
which  are  then  transformed  by  (B4)  and  used  in  (B3)  produce  a 

4*  ~  "  ~ 

There  is,  however,  a  difference  in  the  quality  of  these 

subroutines;  a  survey  paper  by  [26]  may  be  consulted  if  the 
validity  or  efficiency  of  a  subroutine  is  in  question. 
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realization  of  a  homogeneous  impulse  process  with  intensity 
X.  For  our  purposes,  we  set  X°1  and  proceed  to  the  mapping 
step  illustrated  in  Figure  Bl. 

In  his  recent  book41  Snyder  provides  guidance,  but  leaves 
it  as  a  homework  exercise  to  deduce  the  proof  for  a  method  of 
rescaling  the  interarrival  times  of  a  unit  intensity  process 
to  obtain  a  realization  of  a  specified  inhomogeneous  process. 

The  results  are  as  follows:  Let  {t^>  be  the  set  of  occurrence 
times  of  a  unit  intensity  process  as  illustrated  in  Figure  Bl . 
Let  A(t)  be  the  integral  of  the  specified  intensity  function 
X(t) : 

t 

A(t )  -  /  X(a)do  (B5) 

o 

This  function  is  continuous  and  monotonically  increasing  and 
therefore  has  an  inverse  function  A'1  such  that 

-  A^)  (B6) 

tj  -  A*l( tj)  (B7) 

The  set  of  times  <t ± }  generated  by  applying  equation  (B7)  to 
the  set  {t ^ }  is  the  required  realization  of  the  inhomogeneous 
process. 

Algorithm  for  simulation  of  event  times.-  The  subroutines 
which  generate  uniform  random  variables  produce  real  numbers. 
If  the  above  equations  are  applied  exactly,  one  would  be  forced 
to  digitally  solve  an  integral  equation,  thus  producing  another 
set  of  real  numbers  {t^}.  One  could  carry  t»ie  full  resolution 

*See  Snyder  (9]  page  62. 
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of  the  computer  to  the  bitter  ond  of  the  simulation.  However, 
at  some  point  the  use  of  digital  filtering  techniques  to 
simulate  analog  filter  characteristics  of  real  photomultiplier 
tubes  (PMT)  and  other  electronic  devices  would  become  appro¬ 
priate.  At  such  a  point,  the  simulated  signal  would  have  to 
be  Interpolated  and  respecified  on  uniformly  spaced  increments 
with  separation  At.  We  may  greatly  simplify  the  required 
algorithms  by  rounding  the  set  of  occurrence  times  to  the  near¬ 
est  At  interval  and  point.  By  incorporating  this  step  directly 
into  the  solution  of  the  inverse  function  we  avoid 

the  problems  of  solving  the  integral  equation  exactly  and  sim¬ 
plify  that  step  as  well.  The  entire  procedure  is  thus  simpli¬ 
fied  to  the  following. 

1.  Select  a  At  small  enough  to  provide  adequate 
accuracy  for  uniform  sampling  of  X(t)  and  cal¬ 
culating  its  integral  by  the  trapazoidal  rule 
Integration  method. 

2.  Beginning  at  “  tQ  “  0,  compute  realizations 
of  as  discussed  above. 

3.  Calculate  the  trapazoidal  rule  approximation 
of  A(kAt)  for  each  integer  value.  °-5k-kmax* 

4.  Use  conditional  statements  to  test  the  latest 
real  value  of  against  the  integral  A(kAt) 

as  it  is  generated  iteratively  to  determine  a 
histogram  TC  of  the  discretized  occurrence 
times  at  kAt.  It  is  possible,  in  the  simula¬ 
tion,  for  more  than  one  value  of  to  be 
mapped  to  the  discretized  time  kAt. 

There  are  many  ways  that  the  fourth  step  could  be  implemented. 
We  have  elected  the  following:  If  the  value  of  t ^  lies  in  the 

range 

A(kAt)  +  M(k-l)&t]  <  T  <  Aj(k+l)At]  +  A(kAt)  (B8) 

2  — "  1  —  2 

where  A(kAt)  is  the  trapazoidal  approximation  of  the  integral, 
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then  an  occurrence  is  added  to  the  histogram  value  at  k. 

Inclusion  of  random  amplitude  effects.-  When  the  point 
event  represents  a  photoelectron  pulse  or  a  classical  LV  sig¬ 
nal  burst,  we  may  wish  to  assign  the  event  an  amplitude  mark 
(Impulse  weight).  This  may  be  done  by  separately  generating 
realizations  of  additional  random  variables  according  to 
desired  statistics  and  accumulating  one  each  of  these  at  TC 
for  each  occurrence  time.  Note  that  it  is  not  adequate  to 
simply  multiply  each  value  TC  by  a  random  variable,  since  in 
some  realizations  more  than  one  event  contributes  to  the  same 
value  of  k. 

The  statistics  of  the  pulse  height  distribution  of  a  poor 
(PMT)  may  be  nearly  Rayleigh,  while  the  Gaussian  density  with 
15-25%  relative  standard  deviation  may  bo  adequate  to  model  a 
good  PMT.  Very  little  documentation  exists  concerning  the 
probability  density  of  the  classical  signal  bursts.  This  topic 
is  discussed  and  some  data  is  presented  in  our  recent  AEDC 
report  [1] .  Of  the  simplest  densities  an  exponential  density 
or  a  Rayleigh  density  would  be  used  to  simulate  amplitudes  from 
an  unseeded  flow.  As  we  have  shown  it  does  not  follow  that 
the  amplitude  probability  density  agrees  in  any  recognizable  way 
with  the  particle  size  distribution;  even  monosized  particles 
may  produce  a  very  strange  amplitude  probability  density  [1] . 

When  it  is  desirable  to  use  Rayleigh  or  Gaussian  random 
variables  the  following  procedure  is  recommended.  Generation 
of  realizations  of  Rayleigh  or  Gaussian  random  variables  may  be 
obtained  by  first  generating  uniformly  distributed  values  on 
the  interval  (0,1).  Let  x^  and  x2  be  two  such  independent 
realizations.  Then  we  obtain 

R  »  (-2o2  In  xx)1/2  (B9 ) 

8  -  2ttx2  (BIO) 
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where  R  la  a  realization  of  a  Rayleigh  random  variable  with 
parameter  c,  mean  *  ,  and  probability  density  pR(R)  given  as 

pR<R)  “  e”(R2/2a2),  R>0  (Bll) 

'  o 

and  0  is  a  realization  of  a  uniformly  distributed  random  vari¬ 
able  on  the  interval  (0,2tt).  Multiplying  (B9)  by  /2/tt  produces 
Rayleigh  variables  with  mean  "1.  If  Gaussian  random  variables 
are  desired,  the  process  is  continued  from  (B9)  without  the 
/ 2  / tt  factor  by  converting  to  rectangular  co-ordinates: 

X  -  R  cos  6  (B12) 

Y  -  R  sin  0  (B13) 

When  this  is  done,  X  and  Y  are  two  independent  realizations  of 

o 

a  Gaussian  random  variable  with  zero  mean  and  variance  or  ,  i.e., 

Px<x)  “  ~  e“x2/2°2  (B14) 

x  /5? 

and  the  same  form  for  Py(Y).  The  above  procedure  efficiently 
produces  exactly  Gaussian  random  realizations  as  opposed  to  a 
program  such  as  GAUSS  which  sums  12  independent  uniform  random 
variables  to  obtain  approximately  Gaussian  variables  by  the 
central  limit  theorem  (Forrester  [25}). 

Example  Simulation  Program 

We  have  included  at  the  end  of  this  appendix  a  copy  of 
the  printout  of  a  FORTRAN  IV  program  which  is  illustrated  in 
flow  form  in  Figure  B2.  The  occurrence  times  of  the  classical 
signal  bursts  are  generated  as  a  homogeneous  Poisson  process. 

The  amplitudes  of  the  bursts  may  be  either  generated  randomly 
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with  exponential,  Gaussian,  or  Rayleigh  density  with  a  speci¬ 
fied  mean,  or  they  may  be  set  equal  to  a  constant.  The  pro¬ 
gram  computes  the  burst  waveforms  at  each  At  and  sums  all  that 
are  present  to  form  a  classical  signal.  This  signal  is  inte¬ 
grated  in  a  trapazoidal  fashion  and  compared  with  the  occur¬ 
rence  times  of  a  unit  intensity  homogeneous  Poisson  impulse 
process  simulated  as  described  previously.  If  random  amplitudes 
have  been  assigned  to  the  photoelectron  pulses,  these  are  gen¬ 
erated  and  added  to  a  histogram;  otherwise  l's  are  added  to  the 
histogram. 

The  idealized  photon  processor  portion  of  the  program  sums 
the  values  of  the  histogram  (with  no  random  PMT  pulse  height 
effects)  over  an  Interval  At  which  is  some  selected  integral 
number  ( ITAU)  of  At  units  in  length.  The  sum  is  the  photon 
count  sequence  referred  to  in  the  text  as  {n^h  The  ideal 
processor  then  computes  either  the  photon  correlation  for  delay 
values  0-IP  (ONE  ■  false)  or  computes  and  sums  the  dual  corre¬ 
late  and  subtract  terms  (ONE  *  true). 
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APPENDIX  C 


VARIANCE  OF  THE  DUAL  CORRELATE 
AND  SUBTRACT  ESTIMATOR 

In  this  appendix  we  wish  to  determine  a  formula  for  the 
mean-square  deviation  of  the  estimate  ftpq  defined  in  the  text 
equation  (56)  at  value  of  delay  (near  Tm/4  where  the  expected 
value  of  ft__  is  zero.  It  was  our  hope  to  do  this  for  the 

pq 

general  signal  model  which  we  have  presented.  To  complete  such 
a  task  requires  the  use  of  fourth-order  moments  of  the  classi¬ 
cal  signal  process  X(t),  and  it  was  for  this  reason  that  the 
derivation  in  the  last  part  of  Appendix  A  was  undertaken.  We 
did  not  have  time  during  the  contract  period  to  evaluate  and 
use  the  fourth  order  moments  of  the  general  signal.  For  this 
reason  we  have  restricted  this  analysis  to  the  case  where 
steady  background  light  is  the  predominant  source  of  variability 
error.  (Steady  light  adds  variability  error  even  though  it 
cancels  in  the  mean.)  This  simplifies  the  problem  because  of 
the  simplicity  of  the  fourth  order  moment  equations  for  steady 
light  (homogeneous  Poisson  counting  process). 

We  define  ft,  here  as* 

PQ 


ft  ■  1  m, 

pq  1 


(Cl) 


mk  "  nknk+p  “  nknk+q 


where  the  summations  will  all  be  from  1  to  N  unless  otherwise 
noted.  We  have 


<V  '  "V'W  -  <nk>  <nk+q>  <c2> 


*This  derivation  was  performed  with  plus  signs  in  the 
delay  subscripts  instead  of  minus  signs.  There  is  no  differ 
ence  in  the  results. 
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where  n  ■  XAt  ■  nAiP/hv.  Prom  (C2)  we  have 


<fl[  > 

pq 


(C3) 


we  wish  to  evaluate  var  ■  <fi2  >  -  <&_  >^  ■  <ft2  > 

pq  pq  pq  pq 


(C4) 


In  (C4)  the  product  of  sums  was  expanded  and  the  order  of 
expectation  and  summation  interchanged.  The  i  ■  J  terms  are 
separated  because  they  behave  differently. 

The  terms  in  the  first  summation  give 


<m12  > 


-  2<nlnl+pniV  +  <ninl+q>  (CB> 


<n2>l<n2+p>  f  <nf+q>  -  2<i.i+p><B1+q>] 


«  (n  +  n2)(2(n  +  n2)  -  2n2)  -  2n2  +2n3 


where  the  theory  of  homogeneous  Poisson  processes  has  allowed 
the  factorization  due  to  independent  of  nonoverlapping  count 
intervals''1  and  from  which  we  know  that 

<n2>  -  n  +  n2  (C0) 

Next  we  must  evaluate  the  i  /  j  terms  in  (C4).  There 

2 

are  N  -  N  such  terms  but  .nany  of  them  are  zero.  We  have 


* 

2  At  this  point  for  the  general  nonsteady  signal  evaluate 
>  with  X(t)  conditionally  given  and  then  evaluate  fourth 
order  moments  of  the  process  X(t)  which  is  also  Poisson. 
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<iniV 


(C7  ) 


-  ^n^n^  -  <nini+qnjnj+p> 

”  <nini+pnjnj+q>  +  <nini+qnjnJ+q> 

When  all  the  subscripts  are  unequal  we  obtain 

kMA 

<m1mJ>  -iP-rT-rT  +  ri-O 

We  are  restricted  to  1  f  j  and  p  f  q  f  0.  If  we  examine  the 
matrix  of  products  m^nij ,  the  only  allowable  products  for  which 
some  of  the  subscripts  of  n  in  (C7)  are  equal,  excluding  the 
case,  are  found  on  diagonals  parallel  to  the  i-j  diagonal. 
These  diagonals  are  i  *>  j  ±q,  i  ■  j  ±p,  and  i  ■  j  ±  (p-q). 

We  examine  these  diagonals  separately:  for  i  ■  j  +  q 

**1*3*  ‘  <nJ+qnJ+q+PnJn3+P>  •  <nJ+qnJ+2qnjnj+p>  <C8> 

2  ,2 

“  <nd+qnJ+q+pnj>  +  <nd+qnj+2qnj> 

■  n4  -  rT4  -  (n  +  n2)n2  +  (n  +  n2)n2  -  0 

Similarly  for  i  ■  J  -  q,  i  ■  j  ±p,  the  result  is  zero  by  sub¬ 
traction.  For  the  diagonals  i  ■  j  ±  (p-q),  however,  we  obtain 
after  substitution  and  evaluation 

—3 

<mimj>  *  -n  (C0) 

For  large  N  we  may  neglect  end  effects  and  observe  that  there 
are  approximately  2N  terms  which  result  from  these  two  diagonals. 
We  may  now  evaluate  (C4)  using  (C5)  and  (C0)  as 

<Mpq>  “  N<2X“2  +  »3>  “  2Nn3  (CIO) 

-  2Nn2  -  2N(XAt)2  -  Var(ftpq) 

This  is  the  desired  result. 


01 


APPENDIX  D 


SYSTEM  DESIGN  OP  AN  ADVANCED  PHOTON  COUNTING 
PROCESSOR  FOR  LOW-LEVEL  LDV  SIGNALS 

Introduction 

A  problem  confronting  us  in  the  present  design  was  one  of 
complexity  and  cost.  A  two  channel  duplication  of  the  AEDC 
processor  [1]  with  the  extra  requirements  of  dual-channel 
advanced-concept  operation  would  have  been  prohibitively 
expensive.  To  solve  this  problem  we  have  made  several  sig¬ 
nificant  changes,  some  of  which  utilize  the  power  of  a  high¬ 
speed  computer  controller. 

Design  Approach  for  the  Advanced  Processor 

Major  cost  savings  are  associated  with  exact  duplication 
of  circuit  layout.  For  this  reason,  the  dual  channel  system 
is  designed  so  that  it  may  be  operated  either  as  two  com¬ 
pletely  separate  identical  units  in  one  rack;  or  as  two  inde¬ 
pendent  units  with  synchronized  data  (for  cross  correlations) 
or  as  a  one  channel  system  with  a  synchronous  normalizing  chan¬ 
nel  (same  system  clock  and  nk  sequence  but  separate  delays). 

The  counter/timer  (C/T)  functions  needed  for  system  control 
were  obtained  in  the  AEDC  system  [1]  from  a  $1500  laboratory 
counter,  selected  because  of  the  availability  of  options  ($285) 
for  computer  control  and  read.  It  is  considerably  less  expen¬ 
sive  to  include  the  C/T  functions  in  the  special-purpose  hard¬ 
ware  to  avoid  both  the  cost  of  two  units  ($3570)  and  the  asso¬ 
ciated  computer  cables,  I/O  cards,  etc.  This  puts  an  even  more 
stringent  requirement  on  reducing  complexity  of  other  circuits 
to  give  more  room  on  the  wire  wrap  panel  (162-180  IC  sockets 
unless  multiple  panels  are  used).  The  following  system  concepts 
have  been  incorporated. 
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1.  Eliminate  all  front  panel  controls  and  displays. 

The  system  thus  utilizes  the  display  and  command 
of  the  controlling  mini-computer  and  cannot  func¬ 
tion  separately.  This  saves  LED  driver  integrated 
circuits  as  well  as  the  front  panel  itself. 

2.  Use  a  newer,  lower  cost  voltage  variable  oscillator 
for  the  system  clock. 

3.  Replace  the  laboratory  counter  with  ECL  counter/ 
timer  (C/T)  circuits. 

4.  Use  an  external,  separate  heavy  duty  ECL  power  supply 
(“  10  Amps)  to  save  space  in  the  system  enclosure  and 
avoid  a  large  heat  source  in  the  enclosure. 

5.  Reduce  the  maximum  count  to  3  bits  instead  of  4 
(see  following  justification).  This  simplifies 
all  of  the  circuits  and  reduces  package  count. 

6.  Reduce  the  number  of  control  and  read  data  circuits. 
(See  following  discussion. ) 

7.  Do  not  require  the  3  bit  counter  to  be  selectively 
saturable.  Let  it  saturate  at  7  (111). 

8.  Remove  the  single-clipper  circuit  and  the  associated 
multiplexer. 

0.  Limit  the  accumulator  to  15  bits  +  sign  to  be  com¬ 
patible  with  the  18  bit  word  of  the  computer.  Read 
the  accumulator  often  enough  to  avoid  overflow. 

10.  Replace  the  4  bit  subtractor  with  a  3  bit  subtractor 
adder.  The  add  function  could  be  used  in  one  of  the 
normalization  schemes. 

11.  Include  an  experimental  analog  servo  loop  for  zero 
adjusting  the  system  clock  to  the  proper  multiple 
of  the  mean  signal  frequency. 


Reduction  of  Counter  Bit  Number 

The  maximum  periodic  rate  of  discriminator  output  pulses 
is  120  MHz  (dead  time  >  8nsec).  The  largest  random  rate  is 
usually  less,  say  70  MHz.  In  order  for  the  count  of  7  to  be 
executed  therefore,  we  assume  the  count  interval  to  be  larger 
than  0.1  ysec  ( 1 / f Q  >  0.4  or  0.8  ysec  depending  on  the  delay 
choices).  Thus  we  may  have  the  count  exceeded  for  signal  fre¬ 
quency  fQ  <  2.5  MHz  or  1.25  MHz.  We  observe,  however,  that  a 
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count  of  7  in  a  quarter  of  a  cycle  is  a  peak  rate  of  14  elec¬ 
trons/cycle.  This  is  where  the  burst  counter  can  become 
useful.  Thus  a  3  bit  counter  is  adequate  because  when  larger 
signals  are  present  a  burst  counter  processor  should  be  used. 

If  the  larger  signals  are  unwanted  because  they  occur  from 
larger  particles,  then  saturating  the  count  at  7  will  reduce 
their  effect  in  favor  of  smaller  signals. 

Read  logic  for  photon  statistics ♦  -  The  16  bit  word  length 

and  the  high  speed  of  the  H.P.  2100  computer  will  allow  straight 

2 

forward  READ  circuitry.  Each  channel  will  have  a  16  bit  T  L 
latch  connected  to  the  15  bit  plus  sign  accumulator.  In  order 
that  the  counter/timer  and  the  multiple  interval  statistics 
may  be  read,  selectable  steering  gates  will  be  used  so  that  the 
first  6  bits  can  be  connected  to  the  output  of  the  multiplier. 
Sequential  numbers  may  then  be  read  as  follows: 

1.  Set  maximum  delay  in  A+B  and  set  one  side  of  multi¬ 
plier  to  1. 

2.  Run  data  and  clock. 

3.  Stop  clock. 

4.  Read  multiplier  out  with  commanded  single  clock 
advance. 

This  will  produce  a  string  of  more  than  20  sequential  numbers. 

In  order  to  read  the  C/T  total  (8  digits)  32  bits  1b  required 
(BCD  code).  This  will  be  accomplished  by  sequentially  reading 
two  16  bit  words  time-multiplexed  to  the  TTL  output  port. 

Description  of  Circuits 

This  section  describes  the  subsystem  circuits  and  compo¬ 
nents  which  comprise  the  dual  channel  photon  processing  system 
one  channel  of  which  is  shown  in  Figure  Dl. 

External  subsystems.-  The  external  subsystems  are  the 
following : 

1.  Two  photomultiplier  tubes  and  associated  housings 
and  power  supplies. 
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2.  Six  precision  50  0  connecting  cables. 

3.  Two  preamps. 

4.  Two  amplifier /discriminators . 

5.  One  NIM  bin  (rack  mounted). 

6.  Two  voltage  controllable  oscillators  (1-200  MHz) 

7.  Two  rack-mounted  integrated  circuit  power  supplies 
and  connecting  cables. 

8.  A  computer  with  two  microcircuit  I/O  cards  (16  bits 
in  and  16  bits  out  for  each  with  device  command  line 
and  device-ready  flag  line),  and  two  36  pair  twisted 
lead  cables  with  connectors. 

Three  bit  counter.-  This  is  a  dual  section  counter  which 
stops  and  holds  the  count  of  seven  instead  of  turning  over  to 
zero  and  continuing.  Two  counter  sections  alternate  so  that  one 
may  have  data  transferred  and  be  reset  while  the  other  one  is 
counting.  In  order  to  avoid  the  possibility  of  counting  a 
border  line  event  twice,  a  dead  time  between  count  intervals 
approximately  2  nsec  will  be  incorporated.  The  dead  time 
between  input  pulses  (from  the  discriminator)  will  be  .<  10  nsec 
with  a  design  objective  of  8  nsec.  The  alternate  count  inter¬ 
vals  will  be  equal  (design  objective).  A  control  bit  allows 
the  counter  output  to  the  delays  to  be  set  to  zero  (data  gate). 

Delay  sections.-  The  two  delay  sections  A  and  B  are  speci¬ 
fied  as  A  «  (0,1,2,3,4,5,6,8,9,10,11,12}  and  B  -  (0,2, 4, 8). 

This  set  a  .lows  sequential  autocorrelations  with  delays  up  to 
20  in  addition  to  the  dual  correlate  and  subtract  and  the 
normalizing  modes. 

For  dual  channel  operation  with  one  channel  normalizing, 
a  multiplexer  is  provided  which  allows  the  data  from  the  other 
channel  counter  to  be  selected  as  delay  input.  Since  the  sep¬ 
arate  channels  are  identical  physically,  each  has  an  input  and 
an  output  to  the  other  channel. 
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Pour  control  bits  select  the  A  delay.  Two  control  bits 
select  the  B  delay.  Ah  additional  bit  selects  the  output  of 
the  other  data  channel  instead  of  the  counter  output  of  this 
channel. 

Adder /subtractor  &  mode  selector.-  In  the  AEDC  unit  an 
option  for  subtracting  a  constant  from  each  n^  was  Included. 

The  purpose  was  to  ovoid  overflowing  the  accumulator  in  auto¬ 
correlate  mode  when  a  signal  with  large  mean  is  encountered. 

The  ability  of  the  computer  for  high-speed  read  to  DMA 
avoids  the  necessity  of  subtracting  a  constant,  since  the  accu¬ 
mulator  may  be  read  often  co  avoid  overflow.  The  two  control 
bits  allow  selection  of  add,  subtract,  or  add  zero.  The  adder 
subtractor' output  is  four  magnitude  bits  plus  sign. 

Multiplier.-  The  multiplier  has  seven  bits  plus  sign  out. 

A  control  bit  allows  the  undelayed  path  input  to  be  set  to  001 
instead  of  nfc. 

Accumulator . -  The  accumulator  has  15  bits  plus  a  sign 
bit.  It  can  be  reset  by  a  single  pulse.  It  will  be  implemented 
in  2's  complement;  the  computer  software  will  convert  to  sign 
and  magnitude.  The  reset  pulse  must  also  transfer  the  accumu¬ 
lator  values  to  latches  in  the  output  port  to  the  computer. 

Counter/ timer .-  This  section  replaces  the  external  counter/ 
timer  used  in  the  AEDC  System.  It  consists  of  two  8-decade  BCD 
counters  and  associated  input  selectors  and  controls.  Each 
counter  may  select  as  input  either  the  precision  1  MHz  oscil¬ 
lator,  the  system  clock,  or  the  second  discriminator  output 
(external  input),  A  fast  prescaler  selects  divide  by  1,  2,  or 
5  for  counter  #1.  Counter  #2  is  multiplexed-  to  the  16  bit  out¬ 
put  port  (1st  4  significant  decades  and  2nd  4  significant 
decades  separately  selectable. )  The  first  counter  produces 
output  pulses  at  decade  counter  #1  selector.  The  package  count 
does  not  include  multiplexing  and  read  outputs  for  counter  #2. 
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A  gate  generator  block  does  the  following:  It  produces 
the  pulse  which  is  applied  to  the  accumulator  and  the  computer 
flag;  it  programably  reads  and  resets  the  C/T;  and  it  produces 
a  gate  pulse  between  two  successive  counter  #1  output  pulses 
for  application  to  the  counter  #2  control  gates  (counter  #2 
stops  and  holds  when  gate  goes  back  down).  The  programable 
states  include:  a)  read  and  reset  counter  #2  at  the  entry  of 
the  control  word  (CW)  and  at  each  output  pulse  from  counter  #1 
(PFC);  b)  read  and  reset  counter  #2  CW  but  not  PFC;  c)  read 
and  reset  counter  #2  at  PFC  but  not  at  CW;  d)  do  not  read  and 
reset  counter  #2  at  either  PFC  or  CW.  (This  allows,  by  software, 
for  counting  the  photon  rate  over  the  duration  of  an  autocorre¬ 
lation  sequence  and  then  stopping  at  the  end  of  a  number  of 
smaller  intervals.) 

The  control  bits  needed  for  the  C/T  are  as  follows: 

Three  bits  are  used  to  select  which  of  3  inputs  goes  to  each 
of  the  two  decade  counters.  (Assumes  both  will  not  have  same 
input.)  Three  bits  will  select  the  decade  output  of  the  decade 
scaler  (counter  #1).  Two  more  bits  select  the  prescale  divisor 
(1,  2,  5).  Three  bits  have  been  allowed  for  the  four  read/ 
reset  states  (2  bits)  and  an  extra  control  bit  (spare).  Finally, 
there  is  one  bit  which  allows  application  of  the  control  gates 
from  the  other  C/T  to  the  accumulator  read/reset  line  instead 
of  this  channel's  C/T.  (This  is  for  cross  correlation  with 
synchronous  accumulator  read/reset.) 

System  clock.-  The  system  clock  accepts  a  periodic  wave 
form  with  1  positive-going  transition  per  cycle)  and  shapes  this 
into  a  periodic  pulse  train  at  the  ECL  voltage  levels.  It 
includes  buffer  gates  for  proper  fan  out  (5  packages).  It 
Includes  controls  which  allow  the  clock  to  be  stopped  cleanly 
and  a  one-shot  clock  pulse  generator  which  can  be  activated 
by  computer  instruction.  This  feature  allows  sequential  n. 
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values  to  be  read  out  of  the  system  for  multiple  Interval 
statistics.  The  control  bits  are  one  for  on/off  and  one  for 
single  pulse.  There  is  also  an  ungated  clock  line  which  is 
used  as  an  input  to  the  C/T.  There  is  a  control  bit  which 
allows  the  other  system  clock  to  be  selected  (when  the  other 
channel  data  is  used,  for  example). 

Output  port.-  Each  channel  has  one  16  bit  TTL  compatible 
output  port  and  a  1  bit  flag  pulse  line.  The  output  port 
includes  a  multiplexer  (two  16  bit  sections  of  the  C/T,  the  16 
bit  accumulator  output  and  a  7  bit  output  from  the  multiplier), 
a  16  bit  latch,  and  ECL  to  TTL  voltage  level  translators. 

The  output  multiplexer  requires  2  bits  of  control  to  select 
one  of  4  outputs. 

At  the  present  time  it  is  not  clear  which  of  the  following 
approaches  could  be  utilized:  a)  all  ECL  construction  with 
only  17  bits  of  ECL/TTL  translator  as  the  output;  b)  a  separate 
panel  section  of  TTL  circuitry  which  includes  most  of  the  C/T 
and  the  output  latches  for  the  output  port.  The  b)  approach 
would  require  less  expensive  IC's  and  less  power  for  part  of 
the  system.  However,  the  cost  of  a  separate  panel  and  the 
panel  interconnections  may  make  a)  preferable. 

Oscillator  control.-  This  subsystem  provides  automatic 
fine  tuning  Of  the  variable  system  clock  for  the  determination 
of  the  mean  velocity.  The  mean  value  of  the  multiplier  output 
is  negative  while  the  clock  frequency  is  too  high.  Thus  the 
intent  here  is  to  use  fast  digital/analog  conversion  and 
analog  integration  with  controllable  reset,  integration,  and 
hold  states  to  provide  a  control  voltage  to  the  external  volt¬ 
age  controllable  oscillator.  Two  control  bits  are  required. 

This  feature  is  a  research  item.  It  may  later  prove  more 
advisable  to  use  the  D/A  converter  with  a  portion  of  the 
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computer  control  word  bits  so  that  the  mean  velocity  selection 
becomes  a  software  task. 

Input  port.-  All  machine  control  is  accomplished  by  30 
bits  of  information  and  the  time  at  which  certain  of  the  bits 
are  changed.  These  bits  are  latched  into  the  processor  input 
port  15  at  a  time  by  a  single  line  command  pulse  from  the  com¬ 
puter.  All  time-critical  bits  of  control  are  included  in  con¬ 
trol  word  0.  (One  of  the  sixteen  bits  from  the  computer  is 
the  address  of  the  control  word  0  or  1.)  Table  D1  provides  a 
tentative  assignment  of  control  bits.  The  bit  number  refers 
to  the  power  of  2  in  standard  binary  format. 

Package  count . -  An  integrated  circuit  package  count  esti¬ 
mate  of  147  IC's  made  for  one  of  the  two  identical  channels. 

The  estimate  assumes  the  use  of  a  panel  with  162  sockets  and 
design  for  each  subcircuit  with  the  same  approach  previously 
developed  for  the  AEDC  unit.  This  leaves  15  spare  sockets  for 
flexibility  in  design  and/or  additions.  A  slightly  improved 
approach  has  also  been  identified  which  uses  2  bit  arithmetic 
logic  units  (ALU)  instead  of  4  bit  ALU's  to  reduce  circuit 
speed  limitations.  This  approach  would  actually  produce  more 
useable  sockets  because  the  panel  would  have  180  standard 
sockets  instead  of  150  plus  12  4  bit  ALU  sockets. 
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